{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-07-02T18:08:37.746434Z",
     "start_time": "2020-07-02T18:08:36.945475Z"
    }
   },
   "outputs": [],
   "source": [
    "import sys\n",
    "sys.path.append('../')\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "from porousmedialab.batch import Batch\n",
    "import numpy as np\n",
    "import seaborn as sns\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-07-02T18:08:37.750477Z",
     "start_time": "2020-07-02T18:08:37.748407Z"
    }
   },
   "outputs": [],
   "source": [
    "tend = 40 * 24 * 60 * 60\n",
    "dt = 1 * 24 * 60 * 60"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-07-02T18:08:37.756127Z",
     "start_time": "2020-07-02T18:08:37.753600Z"
    }
   },
   "outputs": [],
   "source": [
    "bl = Batch(tend, dt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-07-02T18:08:37.763474Z",
     "start_time": "2020-07-02T18:08:37.758787Z"
    }
   },
   "outputs": [],
   "source": [
    "# ED\n",
    "bl.add_species(name='POC', init_conc=12e-3)\n",
    "bl.add_species(name='CO2', init_conc=2e-3)\n",
    "bl.add_species(name='Fe2', init_conc=0)\n",
    "bl.add_species(name='CH4', init_conc=0)\n",
    "\n",
    "# EA\n",
    "bl.add_species(name='NO3', init_conc=1.5e-3)\n",
    "bl.add_species(name='Fe3', init_conc=17.8e-3)\n",
    "bl.add_species(name='SO4', init_conc=1.7e-3)\n",
    "\n",
    "# Henry law equilibrium:\n",
    "bl.add_species(name='CH4g', init_conc=0)\n",
    "bl.add_partition_equilibrium('CH4', 'CH4g', 1.4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-07-02T18:08:37.768473Z",
     "start_time": "2020-07-02T18:08:37.765213Z"
    }
   },
   "outputs": [],
   "source": [
    "bl.constants['Km_NO3'] = 0.001e-3\n",
    "bl.constants['Km_Fe3_surf'] = 2e-3\n",
    "bl.constants['Km_SO4'] = 0.3e-4\n",
    "bl.constants['k1'] = 1 / 24 / 60 / 60\n",
    "bl.constants['SA'] = 600\n",
    "bl.constants['ro_min'] = 3.84e-6\n",
    "bl.constants['MW_Fe3'] = 106.8\n",
    "bl.constants['Fe3_init'] = 17.8e-3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-07-02T18:08:37.772441Z",
     "start_time": "2020-07-02T18:08:37.769707Z"
    }
   },
   "outputs": [],
   "source": [
    "bl.rates['r_NO3'] = 'k1 * POC * NO3 / (Km_NO3 + NO3)'\n",
    "bl.rates['r_Fe3'] = 'k1 * POC * Fe3 / (Km_Fe3_surf +Fe3) * Km_NO3 / (Km_NO3 + NO3)'\n",
    "bl.rates['r_SO4'] = 'k1 * POC * SO4 / (Km_SO4 + SO4) * Km_Fe3_surf / (Km_Fe3_surf +Fe3) * Km_NO3 / (Km_NO3 + NO3)'\n",
    "bl.rates['r_CH4'] = 'k1 * POC * Km_SO4 / (Km_SO4 + SO4) * Km_Fe3_surf / (Km_Fe3_surf +Fe3) * Km_NO3 / (Km_NO3 + NO3)'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-07-02T18:08:37.777375Z",
     "start_time": "2020-07-02T18:08:37.774068Z"
    }
   },
   "outputs": [],
   "source": [
    "bl.dcdt['POC'] = '- r_NO3 - r_Fe3 - r_SO4 - r_CH4'\n",
    "bl.dcdt['NO3'] = '- 4 / 5 * r_NO3'\n",
    "bl.dcdt['Fe3'] = '- 4 * r_Fe3'\n",
    "bl.dcdt['Fe2'] = '4 * r_Fe3'\n",
    "bl.dcdt['SO4'] = '- 1 / 2 * r_SO4'\n",
    "bl.dcdt['CO2'] = '1 * (r_NO3 + r_Fe3 + r_SO4)'\n",
    "bl.dcdt['CH4'] = '1 / 2 * r_CH4'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-07-02T18:08:37.896743Z",
     "start_time": "2020-07-02T18:08:37.780775Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Simulation started:\n",
      "\t 2020-07-02 20:08:37\n"
     ]
    }
   ],
   "source": [
    "bl.solve()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-07-02T18:08:38.276391Z",
     "start_time": "2020-07-02T18:08:37.898585Z"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "'c' argument looks like a single numeric RGB or RGBA sequence, which should be avoided as value-mapping will have precedence in case its length matches with 'x' & 'y'.  Please use a 2-D array with a single row if you really want to specify the same RGB or RGBA value for all points.\n",
      "'c' argument looks like a single numeric RGB or RGBA sequence, which should be avoided as value-mapping will have precedence in case its length matches with 'x' & 'y'.  Please use a 2-D array with a single row if you really want to specify the same RGB or RGBA value for all points.\n",
      "'c' argument looks like a single numeric RGB or RGBA sequence, which should be avoided as value-mapping will have precedence in case its length matches with 'x' & 'y'.  Please use a 2-D array with a single row if you really want to specify the same RGB or RGBA value for all points.\n",
      "'c' argument looks like a single numeric RGB or RGBA sequence, which should be avoided as value-mapping will have precedence in case its length matches with 'x' & 'y'.  Please use a 2-D array with a single row if you really want to specify the same RGB or RGBA value for all points.\n",
      "'c' argument looks like a single numeric RGB or RGBA sequence, which should be avoided as value-mapping will have precedence in case its length matches with 'x' & 'y'.  Please use a 2-D array with a single row if you really want to specify the same RGB or RGBA value for all points.\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7ff75856e6d0>"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaIAAAD7CAYAAAAo/ZDkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO29eZxUxbn//z69THfPyjADIiLiWhLFDQU0AkLc0LjcRONVLlFzDZh4NTcmMZqoufnGJSYxJsaV3KtGwOVqfqjBuMUFXFHAq6BYgkYUUGDGgdl7Pb8/zume7pmeoWfo6Z6aed6v12ifOnXOeaqbrk8/VU89Zdm2jSAIgiAUC0+xDRAEQRCGNiJEgiAIQlERIRIEQRCKigiRIAiCUFREiARBEISiIkIkCIIgFBXfzioopTzAHcChQBi4SGu9Pu38d4F5QAy4Tmu9RClVCzwAhIDNwIVa61al1A+Bf3Uv/bvW+pdKKQvYCKxzy1/XWl+Vn+YJgiAIudDHvn4scA+OlljAXK21VkpdDvw7sM29fJ7WWnf37J0KEXAmENRaH62UmgLcDJzhGjYKuAw4EggCryilngOuBR7QWt+nlLoSmKeUehyYDUwGbOBlpdRioBVYpbU+LQdbBEEQhP6hL339r4DbtNaPKaVOAm4EvgEcAXxba70ylwfnMjR3LPA0gNb6DdeQJJOAV7XWYa31DmA9cEj6NcBTwPHAZ8DJWuu41joB+IF2YCKwh1LqRaXU35VSKhfDBUEQhLzSl77+R8CTbh0fTp8OTr9+lVLqFaXUTke4cvGIKoEdacdxpZRPax3Lcq4JqOpU3gRUaa2jQJ07FPdb4G2t9Yeu0t6otX5EKXUssBA4qrMRSqm5wFyARYsWTfR4ZHpLEAQhV8LhsH3BBResSiuar7Wen3bc675ea10H4DoQv8PxqgAeAm4HGoHFSqmva62XdGdbLkLUCFSkHXtcw7KdqwC2p5W3pZWhlArijCc2Ad93r1mBM+aI1voVpdQeSilLa52Re8h9w+YDrFq1yj7iiCNyMH3gEYlEACgpKSmyJX3DdPvB/DaYbj+Y3wYT7V+5cmWb1vrIHqr0pa9HKTUDZ25pjjs/ZAF/cD0nlFJPAocD3QpRLm7Fq8Ap7g2nAKvTzr0JTFVKBZVSVcB4YE36NcAsnPkgC3gceEdrPU9rHXfP/wL4T/f+hwKfdhYhQRAEod/pdV/vitAfcaZdVrh1K91z5W6/PxPoca7I2lnS07RIikNwoiIudI1dr7V+wo2kmIsjajdorf+qlNoN+AuOatYB5wEnAg8Cb6Td/irgA5zhuHIcz+gSrfUHPdkkHlHxMN1+ML8NptsP5rfBRPtXrlzZOnHixLLuzvexr38HCABfuLfRWut5Sqk5OMENYeB5rfUverJtp0I0EBEhKh6m2w/mt8F0+8H8Npho/86EqJjkMkdkBNFolI0bN9Le3r7zykUkKfyWZeV8TTAYZMyYMfj9/v4ySxCGFLvaX/Tle1woTOwvBo0Qbdy4kYqKCsaNGzcg/3EkSSQSAOQa9WfbNvX19WzcuJG99967P00ThCHDrvYXvf0eFwpT+4uB9S7uAu3t7dTU1AxoEeoLlmVRU1Mz4D09QTAJ6S8GFoNGiGBgusn5YLC2SxCKyWD9XpnYrkElRIIgCIJ5DJo5omKzfPlyLrnkEv72t7+x++67A/C73/2OffbZh5NOOolbbrmFtWvXYlkWZWVlXHnlley9997E43Guvvpq/vnPf+L1ernxxhsZO3ZskVsjCEJ/s3HjRk4//XQOOuigVNnkyZP5j//4jy51165dy69+9Su8Xi8lJSXcdNNN1NbWFtLcfkWEKI/4/X6uuuoq7r333gz3+JprruHwww/n6quvJpFI8MEHH3DJJZfw8MMPs3z5cgAeeughli9fzo033sidd95ZrCYIglBA9ttvPxYsWLDTetdffz3XXHMN48eP56GHHuLPf/4zV101eDYpGJRCtPil9Tz47Ae0heM7r5wjoYCXc088kH85br9u60yZMoVEIsGiRYv4t3/7NwAaGhr48MMP+f3vf5+qd+CBBzJjxgyeffZZvvnNb3LccccBsHnz5kH1K0cQTGDTY0/w6YMPk8jjBL8nGGTsueewx5mn9/ram2++mbfeegvbtrnggguYNWsWv//97xk5ciQA8XicQCCQN1sHAoNSiB5buj6vIgTQFo7z2NL1PQoRwH/9139x9tlnc+yxxwJOmOeee+7Zpd6ee+7J5s2bAfD5fPz0pz/lueee49Zbb82r3YIg9Mymx57IqwgBJNrb2fTYEzsVovXr1zNnzpzU8dlnn83GjRt56KGHCIfDfOtb3+KrX/1qSoRWrVrFwoULWbRoUV7tLTaDMljhzOn7EQp483rPUMDLmdN7FiGA6upqfvazn3HllVeSSCSIRqMpwUlnw4YNqbkkgJtuuolnnnmGa665htbW1rzaLghC9+xx5ul4gsG83tMTDObkDSWH5pJ/W7Zs4b333mPOnDlcdNFFxGKxVP/x97//nV/84hfMnz+f4cOH59XeYjMoPaJ/OW6/nXou/cnMmTN57rnnWLx4MT/5yU8YO3YsixYtYvbs2QC89957vPDCC3zve9/jscceY8uWLcybN49QKIRlWXi9+RVRQRC6Z48zT+/1EFp/LWjdZ599mDx5Mr/61a9IJBLccccdjBkzhscff5yHH36YBQsWMGzYsLw+cyAwKIVoIPDzn/+cN95w8rvedNNN/OY3v+Hss8/G4/FQWVnJHXfcQWVlJSeeeCJXXXUVs2fPJhaL8bOf/WzQjf8KgpAbM2fO5M033+S8886jtbWV448/nlAoxPXXX8/uu+/OpZdeCsBRRx3FZZddVmRr88egSXq6du1axo8fXySLcqevv6QGSvtMTPbYGdPbYLr9UPw27Or3aaCm+EmSrX0DOenpwHwXBUEQhCGDCJEgCIJQVESIBEEQhKIiQiQIgiAUFREiQRAEoaiIEAmCIAhFRdYR5ZH58+fz2muv4fF4sCyLH/7whxx88ME89dRTLFy4EI/HQywW4+yzz+Yb3/hGxrXXXHMNVVVV/PjHPy6S9YIgFJJ169bx29/+lra2NlpbW5k+fTqXXnopDQ0N3HTTTWzevJl4PM7uu+/OlVdeyYgRI2hqauInP/kJzc3NRKNRrrzySg4//PBiN2WXESHKE+vXr+eFF17gwQcfxLIs1q5dy09/+lOuuOIKHnroIe666y4qKipobW3lBz/4AaFQiFmzZgFO5u0PP/yQo446qsitEAShEDQ2NnL55Zfzpz/9iXHjxhGPx/nBD37Agw8+yJIlS/jOd77D8ccfD8Brr73GvHnzeOSRR7j33nuZMmUKF1xwAR9//DE/+tGPWLx4cZFbs+sMSiHa/sYTNLz8MHYkf4kMrZIg1VPPYdiU7KlAhg8fzubNm3n00UeZNm0a48eP59FHH+XSSy/lxz/+MRUVFQAEg0GuuOIKfvnLXzJr1izefvtt3nnnHc455xw+/vjjvNkrCEJuvP7SRyx99kMieUyUXBLwMv3EAzj6uH2znn/++eeZPHky48aNA8Dr9XLTTTfx0UcfsXTp0pQIARxzzDGMHTuWt956iwsuuCC1CHgwZeEelHNEO5Y/kVcRArAj7exY/kS354cPH86dd97JqlWrOOecczj55JN58cUX+eyzz7psdJfMvL1161Zuu+02rr322rzaKghC7ry+9OO8ihBAJBzn9aXd/7DcunVrl6z8ZWVlbNy4scds/ZWVlQSDQbZt28ZPfvITLr/88rzaXSwGpUdUNfn0fvGIqiZ3nxhxw4YNlJeXc+ONNwKwevVq5s6dy4EHHsimTZuoqqpK1f3kk0/Yfffdefrpp2loaGDu3Lls27aN9vZ29tlnny7zR4Ig9B9HT9+nXzyio6fv0+350aNH8/7772eUffbZZ9TW1rJp06Yu9Tds2MAxxxwDgNaayy+/nCuuuIJJkyblzeaiYtu2cX8rV660O/P+++93KSskzzzzjH3BBRfY7e3ttm3bdlNTk33CCSfYzz33nH3++efbTU1Ntm3bdmNjo/3v//7v9pNPPplx/V//+lf7t7/9bbf3L3b7koTDYTscDhfbjF3C9DaYbr9tF78Nu/p9isfjdjwe7/P1TU1N9qmnnmpv2LDBtm3bjkQi9ve//3174cKF9llnnWU///zzqbpLly61Tz/9dDsWi9nr1q2zTzrpJHvt2rU93j9b+1asWNFiD4D+O9vfoPSIisGJJ57IRx99xNlnn01paSm2bXPFFVdw/PHH09raykUXXYRlWcTjcc466yxOOeWUYpssCEKRKC8v59e//jVXX301tm3T0tLCjBkzOO+88zj55JO54YYbuPvuuwEYNWoU8+fPx+v1cvPNNxOJRLj++utT97nzzjuL2ZS8INm3C4xk3y4+prfBdPuh+G2Q7NsDi4H5LgqCIAhDBhEiQRAEoaiIEAmCMCQxcVoiF0xslwiRIAhDjmAwSH19vZGddk/Ytk19fT3BYLDYpvQKiZoTBGHIMWbMGDZu3Mi2bdv6dH1SwCzLyqdZeSEYDDJmzJhim9ErRIgEQRhy+P1+9t577z5fX+yov8GGDM3lkXXr1jF37lzmzJnDN7/5TW699VY+++wzvvWtb2XUe+ihh/jTn/6UOk4kElx00UU8+OCDhTZZEASh6IhHlCe6y6b7yiuv7PTaP/zhD+zYsaMAVgqCIAw8BqUQ/e2Df/DIe0toj4Xzds+gL8DZB32d0w48Puv57rLpbt26tcc07U8//TSWZTFt2rS82SoIgmASOxUipZQHuAM4FAgDF2mt16ed/y4wD4gB12mtlyilaoEHgBCwGbhQa92qlPoh8K/upX/XWv9SKRUCFgIjgSbgfK1132YQXZbof+RVhADaY2GW6H90K0TdZdP1+/2sX7+eOXPmAM4k59atWznttNP48MMPWbJkCbfeeiu33357Xu0VBEHoDX3s68cC9+BoiQXM1VprpdRpwLVu3Xu01n/u6dm5zBGdCQS11kcDVwI3pxk2CrgM+CpwEnCjUirgGvCA1noq8DYwTym1DzAbOAY4GjhRKXUI8D1gtVv3fuDqHGzqka+r4wn68rtPR9AX4OsquwiBk033iy++yCj77LPP+Pzzz9lvv/1YsGABCxYs4P777+eCCy4A4LHHHmPLli2cf/75LF68mPvuu49ly5bl1W5BEIQc6Utf/yvgNq31ccANbrkfuAU4EZgOzHWv75ZchuaOBZ4G0Fq/oZQ6Mu3cJOBVrXUYCCul1gOHuNfc4NZ5yn19G3Cy1jruNswPtLt1f5NW95ocbEpFrSSxbTuV/+nUA2Zy6gEzc7lNr0k+ozPTp0/nrrvu4pxzzmHs2LFEo1FuvPFGjjnmmAzb0jPOpm8Lftttt1FbW8uxxx6b9Rm2bXdpczGIRqPFNmGXMb0NptsP5rfBdPu7oS99/Y+A5AS3D6dPHw+s11o3ACilXgGmAo909+BchKgy7UEAcaWUT2sdy3KuCajqVN4EVGmto0CdUsoCfgu8rbX+UCnVpW42I5RSc4G5AAsXLszB7MKSzKZ77bXXkkgkUtl0p02bxuOPP15s8wRBGOI0Njb6lFIr0orma63npx33uq/XWtcBKKUU8Dscr2pEtro92ZaLEDUCFWnHHtewbOcqgO1p5W1pZSilgjjjiU3A97PcI1W3M+4bNh+c7Nud4/ctyyp6JtwJEyZw//33dyl/5JGOHwKJRIJzzz23i62XXXZZj/e2LGtArVkYSLb0FRPbsHXpMjbcv4hIfR2B2lrGzpnNyOnmBrqY+BmkY5L9lZWVMa31kT1U6Utfj1JqBs7c0hx3fijQXd3uyKXnfhU4xX3gFGB12rk3galKqaBSqgrHJVuTfg0wC3jZ9YQeB97RWs9LDtFlq5uDTYIw5Ni6dBkf3X4Xkbo6sCG8rY6Pbr+LrUtlXlHIC73u610R+iPOtEvS21oL7K+UGq6UKgGmAa/39OBcPKLFwAlKqddwoiIuVEpdjjMG+IRS6lYc8fAAP9datyulrgP+4kZZ1AHn4bhs04GAUmqWe++rgDvduq8AEbeuIAid+HTBIhLhzGjQRDjMpwsWGe0VDWa2Ll3GpwsWEa6rJ1BbM9A92L709X8ASnD6cACttZ7nXveMW/cerXXX/c/TGFQb4x144IEDMvdTOn3ZUMu2bT744APZGC9PmNqGV888C7J9Xy2Lrz72aOEN2gVM/QyS5GJ/0oNN//HgCQTY95KLiyJGsjFeAZBsusJgJ1Bb06vygcrWpct453uX8ta3zmXFRfMG7dBiTx6skMmgyaywq9l0C0VfsvaamE1XyD9j58zO+gt77JzZRbSqd3T2EpLzXMBAHrLqE+G6+l6VD2UGjRDtajbdQmH6kEQSw8a++8RAa2Py2U7U3MCwqbcMpXmuQG0N4W11WcuFTAaNEAkd9HcHWvfyK2y4+78H9a/agfrLfeT0aQw7egpg5o+ZoeQlDAYPtlAMmjkiwSHZgYa31YFt90uI76YHHh70Y98yvt8/DJZ5rlwYOX0a+15yMYERtWBZBEbUFi1QYaAjHtEgoxBDH5H6rsMNMLh+1Q6lX+6FZKh5CSOnTxPhyQHxiAYZhehAS2pqs5YPpl+1Q+mXeyFJegklteIlCB2IRzTIKMQE6R7nnZMxRwSD71ftUPvlXkhMn+cS8o94RIOMsXNm4wlkboGR7w60duqxg37sW8b3BaFwiEc0yEh2lP0ddjwUxr6HQhsFYSAgQjQIkQ5UEASTkKE5QRAEoagMOY+odeMmGt5aQSIWcxJIuruldnndT8Tjzu4XXq+3357Rn5huP5jfBtPtB/PbYKT941WxLeiWISVE25a9zLo//Ak7Ht95ZUEQhEFE8NqfFduEbhkyQvT5U0/z8d3/3a/ejiAIgtB7Br0Q2bbNxkf+yqeLHkyVhcaMoWbKJOfAssCynGzYnV/3A/G4s/Ou12vmW2+6/WB+G0y3H8xvg4n2by22AT1gzrvYB2zb5pP77mfzY0+kysr335+vXPtz/JUVPVzZf5iefdt0+8H8NphuP5jfBhPt37pyZbFN6BYjo+bC0Tjfue5ZXlr5Wbd17Hic9bfdkSFCVYdM4KD/94uiiZAgCILQFWM9om0Nbdz2yDsAHDdxz4xziWiUD2++hfrXl6fKhk+ehPrxD/EY9AtGEARhKGCkEFk4AQfhaJz7n1qbIUSx1jY+uPEmdry7OlU2cuYM9vuP72GZFGopDLiN6QRB6B+MFCI/CbzEieOlrqEtVR5va+O9a39J87p1qbLRp3+dcReej+XJbRSyac0yGl5cRKyxHl9lDdUzZlNxsHR+hWagbkwnCEL+MXKOCGxqPU0A1FaHUqVb/vF8hgiNnX0u475zQa9EqO7Ju4g11gE2scY66p68i6Y1+dtUTsgN2ZhOEIYOhgoR1HibCfi9fHvW+FRZ62cbU6/3+Jcz2PNbZzmh2DnS8OIi7Fhm52fHwjS8KJ1foZGN6QRh6GDk0BzAnPJXaD10NJPS5ofaN29IvW5b9zxNa8b1algt1pi9k+uuvD+QeRGHQuyrJAjCwMBYj6jUijDi/UdSw2ZNa5bR/tn61Hk73NzrYTVfZfZOrrvyfJOcFwlvqwPbTs2LbF069IYGC7GvkiAIAwNjhQgyh80aXlxEIpZInfP4ej+sVj1jNpYvs/OzfAGqZxSm85N5kQ5kYzpBGDoYOzSXJDlsFmusx07LZWp5M8/nQnIYr1hRczIvkonsqyQIQwPjhSg5bOarrCER75hTSApRb4fVKg6eVrRwbZkXEQRhKGL00BzektSwWfWM2Rkekcdb2GG1fCDzIoIgDEWMFqJhx5yZ8l7Kv3IsdscUEf5hNdSeerFRi1FlXkQQhKGI0UNz/uG7p17H29txM//gCQbZ6wfzu9Q3IWuCzIsIgjDUMFqIYg1bUq/jrR2pfryhUJe6yawJyQWryawJwIATI0EQhKGE0UNz0e3pQtSSeu0rK+1SV7ImCIIgDEzM9ojShCiW7hGVdhWiQmVNkMwIgiAIvcNsjyhjaK419dqXRYgKkTVBMiMIgiD0np16REopD3AHcCgQBi7SWq9PO/9dYB4QA67TWi9RStUCDwAhYDNwoda61a0/AngNmKC1bldKWcBGIJk2+3Wt9VW5GB9v+pJELILHV0KspUOIsnlE1TNmZ8wRQf7Du3vKjCBekSAIA5m+9PVp5/4TGKW1vtI9vhz4d2CbW2We1lp39+xchubOBIJa66OVUlOAm4Ez3IeNAi4DjgSCwCtKqeeAa4EHtNb3KaWudI2/RSl1EvBrYLe0++8LrNJan5aDLQ6pjNo2sR3bKKnZg3hbuhB1DVYoRNYEyYwgCILB9KWv9wB/BiYDf0271xHAt7XWK3N5cC5CdCzwNIDW+g2l1JFp5yYBr2qtw0BYKbUeOMS95ga3zlPu61uABHA8kG7cRGAPpdSLQBvww56UE8CmY2uHtm2boGIE4camVJkVCBKJRLpcFzhgCqMOmJJRlq1eXympqSFS1zUzQklNTeo50Wg0b88rBqbbD+a3wXT7wfw2mG5/N/Slr18P3A/8Azgwrf5E4CpXwJ7UWt/Y04NzmSOqBHakHceVUr5uzjUBVZ3Kk2VorZ/TWnd2Dz4HbtRaz8ARrIXZjFBKzVVKrVBKrbDtNGO2bwUgkRGs0NUjKgS10w7H6vSOWh6nXBAEoZg0Njb6kn2o+ze3U5Ve9/Va6wat9bNZHvcQcDEwEzhWKfX1nmzLxSNqBCrSjj1a61g35yqA7WnlbWll3bECZ8wRrfUrSqk9lFKW1tpOr6S1ng/MB3hz+fLUObu5jpKSEuxwe6puSUUFJSUlOTQtv3gb36ZyHDRthEQEPCVQMcYp72xPMezLJ6bbD+a3wXT7wfw2mGR/ZWVlTGt9ZA9V+tLXd8Gd9/+D1nqHe/wkcDiwJFt9yE2IXgVOA/7XHTdcnXbuTeB6pVQQCADjgTXuNacA9wGzgJd7uP8vgHrgN0qpQ4FPO4tQZ9KH5qKuR5QRNZdlHVEhiDXWU1oLpbVdywVBEAY4fenrs1EJrFFKjQdacLyie3p6cC5CtBg4QSn1GmABF7oREeu11k8opW7FERoP8HM3Eu464C9ulEUdcF4P9/81sFApdSqOZ3TBzgxKV6lkdoWMdUSh4giRr7KGWGPXOaJCbawnCIKwC/S6r892E631DqXUz4AXcaLvntda/72nB1u23aPzMSB5Y/lb9sh//BoAqyTEuB8v4L1rf8mOdx0BP+iX1zLssEMLblfnNELghIinJ19NBi2Y5NKnY7r9YH4bTLcfzG+DifavXLmydeLEiWXFtiMbRmZWsIG4N4A3HsaOtJFoa8oYmsuWa64QFHtjPUEQBBMxUogAwoHhlLZ+DjgZFmLpQlSkOSIo7sZ6giAIJmJsip9W37DU69j2LRnZt32lA9L7FARBELJgrBA1eqtSr6Pbt2QOzRVpHZEgCILQe4wVou12R0h7tO5zEskMCR5Pl+22BUEQhIGLsXNEdYkOIQrXfZ567SstxbKsbJcYgQm7yAqCIOQTY4VoS7RjHihavzX1OlvmbVOQXWQFQRiKGDs090V7ENwMC9EdDalyk+eHZBdZQRCGIsYK0Y42G2/FcAASsY5Fudk2xTOFQu0iKwiCMJAwVoha2qP4hznbGtnxjvJiriHaVQqxi6wgCMJAw1ghSiRsrMoRzut0ISpSnrl8UD1jNpYvM+Iv37vICoIgDDSMDVYAiJc5aa7TPaJiZd7OB5IiSBCEoYjRQhQNOUNWGUNzRcozly8kRZAgCEMNY4fmANr8TpqfRIZHJOl9BEEQTMJoIWr2OkI0mDwiQRCEoYbRQtSUCGL5AyRiHWUmR80JgiAMRYwWopb2GL5hu2UGKxi8jkgQBGEoYrYQtTlriTLCt0WIBEEQjMJoIWpui+KrzvSIRIgEQRDMwngh8ncZmpNgBUEQBJMwWoiyD81J+LYgCIJJGC1EzW1RvMNGdhqaE49IEATBJIzOrNDSFsUTqOwo8IDl6XlTvPff2czad78gkUgAYNs9Vs87yed6PGb+BjDdfjC/DabbD+a3wUT7950wcDcMNVqImtuiJKIdi4g8XmczuWRW7s5s2dzIowtWQYHFRxAEodjsO2H3YpvQLebIeRZa2qLEW9tSx5YXYg1buq3/xrKPRYQEQRAGGEZ7RG3hGJGm5tSxxwvR7VvINkvU3NjOmlWbU8cnnXkQFZXOlguWVTiXNRZzPDifr//e+rZP36dlzTLibU14QxWUHTyN0Niv5OXehbC/vzG1DU8+uprWlkiX8tKyEk49a0IRLOo7pn4GSUy0vzW6eeeVioQ572Ia6bLRsqNDiCwvxLZn94jeem0D8bgzrjtmr2omT927P03slkjE6UhKSkr65f5Na5ZR997/MCwehhIgDtZ766gdd3Fesnr3t/2FwNQ2PPKXlVnLW1sijD9k4A67ZMPUzyCJifavXClClGc6pKh1e2PqtccL0SxDc9FonJWvbUgdT5leHBEqBA0vLsKOhTPK7FiYhhcXGbm9xOqVG3nhKc2OhjaqqkPMnKWYMHHMkHl+OlXVIXY0tGUtFwSTMXOOKM0lamts6ijuxiNavXJjakijqjrEgQeP6ncTi0Wssb5X5QOZ1Ss3suSR1anOd0dDG0seWc3qlRuHxPM7M3OWwu/3ZpT5/V5mzlJFsUcQ8oWRHlH60Fx7Y0tHuTtHlI5t27yx7J+p40lT98bjNVN/c8FXWUOssS5ruWm88JQmGo1nlEWjcV54ShfEKyn28zuTfObzf/+Axu3tRffQhJ0zkDzqgYyRQpROpLklJUweLyTamom3t+ANOhkWPtLbqNvizCOVBLwcPmnPIllaGKpnzKbuybsyhucsX4DqGbOLaFXfyDYM1VP5YHt+NiZMHIOaMBIwa35iKJL0qJM/ZpIeNSBi1AkjXYP0KLdoS2tHuTtqkT4898bSDm/o8EljCYb8/W9gEak4eBq1p16Mr7IWsPBV1lJ7an4CFQpNd3MfhZoTKfbzBbPpyaMWMjHeI4q1tJD8XehxhSi6fQuBUfuw9fNGPv5wGwCWBZOmjv33CaUAACAASURBVCuKjYWm4uBpRgpPZ2bOUhm/KKGwcyLFfv5gZfXKjUNieHEgetQDFSOFKH2OKNGWuaAVOha1Lk+bG1IHj6K6RhKimkSycyrWGHuxnz8YGUrDVRLlmDtGClG6EtlpQpTuEbU0hXl31abUuSnT9ymUdUIemTBxTFE7qGI/f7Ax0AJA+hPxqHPHSCHKyIMQbu8od1sT276FFa9tIB5zFrCO3rOKPcdV53TvrUuX8emCRYTr6gnU1jB2zmxGTjd/mEsQBgJDabhKPOrc2akQKaU8wB3AoUAYuEhrvT7t/HeBeUAMuE5rvUQpVQs8AISAzcCFWutWt/4I4DVggta6XSkVAhYCI4Em4Hyt9bYejUpTIk8kLTrM9Yja6rex4oNPUuVTpu2TUxqfrUuX8dHtd5EIO/cMb6vjo9vvAhAxEoQ8MNSGq0zyqPvS16ed+09glNb6Svf4NOBat+49Wus/9/TsXKLmzgSCWuujgSuBm9MePgq4DPgqcBJwo1Iq4BrwgNZ6KvC2azxKqZOAZ4H09NjfA1a7de8Hrt6ZQVaaEnmjHR5Rcmjuo23ltDQ7C1grq4KMPzS39CefLliUEqEkiXCYTxcsyul6QRB6RhblDmh63dcrpUJKqYXAJWl1/cAtwInAdGCue3235DI0dyzwNIDW+g2l1JFp5yYBr2qtw0BYKbUeOMS95ga3zlPu61uABHA8kJ4061jgN2l1r8nBphTeWEcSSG95JXa4kbVt41NlRxy9J/F4jHg829WZhOu6LgR1yutTuaV2lWg0mpf7FAvT7Qfz22Cy/WrCSGKx8bz0zDqadoSpHBZk+on7oSaMzNt3rBCY/Bn0QF/6+vU4DsQ/gAPduuOB9VrrBgCl1CvAVOCR7h6cixBVAjvSjuNKKZ/WOpblXBNQ1ak8WYbW+jnXsO7un6rbGaXUXGAuwIIFCwHw2An8CXc/IsvCXzOSTz8tY3t8OAB+v4fDJuXuFpfU1BLJIkYlNeZlJehv3nt7M0ufXU/j9vZUZ3LQ4aOLbZZgAAcdPpoDDh4BgN8/uNf1DSQaGxt9SqkVaUXztdbz04573de7YvOsUuqCHu7Tbb+eJBchagQq0o49rmHZzlUA29PK29LKcrl/t3XdN2w+wKpVq2yAkkSaN1QaIlAzmrUfBlJlh00aS2VV7iHbe317dsYcEYAnEGCvb8/O+yp2k1fFv/f2Zp5evDYVDdS4vZ2nF6/F5/MZMx4OZn8GYL79YH4bTLK/srIyprU+socqfenrc7nPzjQgpzmiV4FTAJRSU4DVaefeBKYqpYJKqSocl2xN+jXALODlXO6fQ90Ufp+HQKLDPfaVltJojWRzNNkR2kye1rss2yOnT2PfSy4mMKIWLIvAiFr2veRiCVToxNJn18uKcUEYfPSlr8/GWmB/pdRwpVQJMA14vacH5+IRLQZOUEq9hhOvdqFS6nKcMcAnlFK34oiHB/i5Gwl3HfAXN8qiDjivh/vf6dZ9BYjspG6KspCfkpYOIfKWlrJuW4cIj6ttZ3ht7xewjpw+TYRnJzRub89aPhhDcAVhCNHrvj7bTbTWUfe6Z9y692itN2Wrm8SybfP2zl61apX95+casDZ8xOxNzwBQMf5A3hsznXX/dIbrjtvnM6Zd8r1impkVEzfUSicSiXDHTcuyilFVdYgfXP21IljVOwbDZwDm2g/mt8FE+1euXNk6ceLEAZlexsgFreB4RPG0OSJfaSmNLR1h3eW2efvv5Iv+Tj0//cT9MuaIQEJwBUHoO8YKUXnIT1sifWguROP2jo4xGMm+ZfhgpxC5vA46fDQ+n09WjAuCkBeMFaKykJ9EmkeUCJbT1uYEeHiIE2j7AjsRx/J4u7vFoKRQubxMWjEuCMLAxsj9iMARokC8wyNq95V3nPM0Y5Eg3tqU7dJBzVDK5SUIwuDAWCEqD/kzwrdbrWDHOY+zI2u8uaHgdhUb2cxNEATTMFyIOobmWu2OhazlXscTijd/WXC7io3k8hIEwTQMniMqoSTNI2qO+QDnuMz1iGIGekS7GvEmqecFQTANY4WoPOQnmC5E4bTQ7dTQXI9ZJQYc+Yp4k0ACQRBMYtAMzTW1dSzM7RiaM8sj6iniTRAEYbBirBCVhfypoTkbaGzp6MCTHlGsyaw5Iol4EwRhKGK0ECWj5mKeEiIRZ1twv88iYDnpZ0zziCTiTRCEoYixQlRe2jE01+bvSHZaVR0guSu4aUIkEW+CIAxFjA1WCAV8KY8ofTHrsJpy2Oq8jjVvx7ZtLMvKdosBh0S8CYIwFDFWiKxYFA9OgEKLvzJVXl1TjrWjFDvcCokYibYmvKWV3d1mwGFKxFvTmmU0vLiIWGM9vsoaqmfMpuJg2T5DEITeY+zQXLylNfW6paRDaIYND+Err+6oZ9jwnAm0vP8KdU/eRayxDrCJNdZR9+RdNK1ZVmzTBEEwEGM9olhrS+p1m69jjmjY8FK85dVE6519mGLNDZSM3KugtvX3NgzFpvHlh7Fj4YwyOxam4cVF4hUJgtBrjBWieGtHSHO7P22OyBWiVL0Ch3AXYhuGYhNvrMtaHmscuntACYLQd8wdmmt1huZsIOotTZUPGx7KEKJYgbMrDIVFqd7K2qzlvsqaAlsiCMJgwFghirlCFPEGsd09hwJBH6HSErwVxZsjGgqLUiunnoPlC2SUWb4A1TNmF8kiQRBMxuChOUeI2tPmh6qHO55RZrBCYYfmqqpDWUVnMC1KLfvKsfh8PomaEwQhLxgsRE5n35Y2P1Q13OnsM4fmCusRzZylMuaIYHAuSq04eJoIjyAIecFYIUoOzXUOVIBMISr00JwsShUEQegdxgpRvMUJ3+4cug3gKx/eUa8I2RV2ZVFqwk4QjkVojbYRjoWJJeLE7QSxRIx4Iu4eu/9PxLGxsW2bhG1jk+h4bdupc86yX/e/XY47np08A2mFnYjFHE/P5+tIRWR3X31AEovHAPB5zfznb7r9YH4bTLS/lvKdVyoS5ryLnYi5Q3PZPCJPIITlD2JH27FjERLtLXhDxf0QIrEInzdv5dOGTWxu2sq21nqaIs20RtpojbXTGm2jLdpOW7Q9TRAEQRDyw0/3u6jYJnSLsUKUDFZoywhW6AgI8JYPI9bwhVO3uaGgQrS1pZ7/+/w9Njd+waamLWxu2kJdy5ciMIIgCFkwWohsoN1flipLekQAvorhGULEiD371Z6EnWD1lg94et1LrNq8ZpdEJ+ALUOoLEvCV4PP48Hq8+CwvPo8Xr/vn8/jwWh48lgfLsrAsCw+W+9qDBwss8CQj9C2wsJIvwbIyj5NYWcrS25lwttvweDIj/61urxh4xBPO8KLX491JzYGJ6faD+W3Ixf5lG5bT3ikDCUDQF2DaXpP7zTYTMVaIYq2thL2l2JbzDyFuQUmgozmZkXP9F8LdGmnjpU9e55n1S/m8aWu39SzLYreyWkaVjWB0xW6MGTaa4aEqSv0hQv4gIX+IUn+QkC84oL+ckYiz9UZJSUmRLek7prfBdPvB/DbkYv9zH2XPvRiOhbnoyHP7xa6eWLlyZcGfmSvGClG8tTVjH6LOvzsyI+fyn13hsx2beWbdUpZuWE44y6+eQ0eN5ysjDmB05W7sUTGK3cpr8Xv9xn8BBUHIjZrS4dS1dv0RXFM6PEvtoY3BQtRGu79DbNrtBNFYHL8bzeXrp7VEjeFmbl9+H29//l6XcyF/kBl7H8OJ+01jdMVueXumIAjmce4hZ3D3W4uIxCOpshJvCececkYRrRqYGCtEsdZW2kId8z5hoLktSnWFI0T9sZYoEovwm5fv5MP6jzPK96wazcn7HcfUvY4i6A/m5VmCIJjN1L0mAfDgu49T3/olNaXDOfeQM1LlQgfGClGivZ32yo5IuDDQ3BqlusIRAl+eM3An7AR/Wn5fSoQsLCaPOZyT95/O+BH7G7MLrCAIhWPqXpNEeHLATCFyV1Cmh26HsWlpj6aO853mZ+E7i1m+8e3U8fmHn8UpB8zc5fsKgiAMdczMvu2GEKcvZo3geERJOgcr2Luw/P/pdS+xRP8jdXzK/jNEhARBEPKEsUKUwCLs61hDFAZa2jqEyBMsw/I5kWl2tB070rdtGFZseod73/7f1PGkPQ7j24ed1Te7BUEQhC4YOTRn2zZhXxm25ehoxF0+mj40Z1mWk11hu7O2J9bcwPY3VvDpgkWE6+oJ1NYwds5sRk7vPoP0+vpP+MPr/5PypvYfPo5Lp1zYZTGnIAiC0Hd2KkRKKQ9wB3AojuNxkdZ6fdr57wLzgBhwndZ6iVKqFngACAGbgQu11q3d1B0OfAiscW+5WGv9xx6NSiSyriFKH5oDZ3guKURbX3yJz/73SRJhp3Z4Wx0f3X4XQIYYrV65MZU5OxZoJ7THCCK1m9mtrJYrpn6PgE/W/wiCMPjIc19/K/BVoMm9/Ayt9Y7unp3LT/szgaDW+mjgSuDmNMNGAZe5DzwJuFEpFQCuBR7QWk8F3gbm9VD3COBBrfVx7l/PIgRg27T70ueHHI8lfWgOMrNwb/7bcykRSpIIh/l0waLU8eqVG1nyyOrUxna+cJA9PpnAbtv35qrp/0FVsHKnpgmCIBhKXvp695IjgJPS+vVuRQhyG5o7FngaQGv9hlLqyLRzk4BXtdZhIKyUWg8c4l5zg1vnKff1R93UnQgcoZRaCmwFLtNaf96jRYlExoZ4SXlpbGlPZS4AsEo7hCO6vYlshOvqU9c8//cPMja0A/AkfIz9/BBqA9UZ9+4r0Wh055UGMKbbD+a3wXT7wfw2mG5/N+Slr1dK/RHYH5ivlNoN+B+t9T09PTgXj6gSSFezuFLK1825JqCqU3m2svTyD4BfaK2nA48Bf8pmhFJqrlJqhVJqRSIez/CIkkLU0hbLbFzZsNRrX3n2rbpLampSrxu3t2et09YUy1ouCIJgCo2Njb5kH+r+ze1UJV99fRlOP/5vwMnA95VSh/RkWy4eUSNQkXbs0VrHujlXAWxPK2/LUta57nKg1S1bDPy/bEZorecD8wFWvvSSnTlH5AzNtYZjGTncAlW1qdc1R4xl2/JPMobnPIEAe317duqaymHBrGJUVR3Ke24403PNmW4/mN8G0+0H89tgkv2VlZUxrfWRPVTJV1/fCvxRa90KoJR6AWfe6d3uHpyLR/QqcIp7wynA6rRzbwJTlVJBpVQVMB4n6CB1DTALeLmHuv8NfNOt+zVgpyli7YTdZQ0RdJ0j8lZ0zBGV7uZj30suJjCiFiyLwIha9r3k4oxAhX2OKSPhyfR+/H4vM2epnZkkCIJgOvnq6w8AXlFKeZVSfpzhu1U9PTgXj2gxcIJS6jWcbWouVEpdDqzXWj/hRke8jCNqP9datyulrgP+4kZZ1AHnaa1buql7JXCPUur7QAuw020Ek+HbuAZ1CFGmiPg65ZsbPX1aj+HamyrWs2ncJ+y2UVESCVFVXcrMWarP234LgiAYRD77+kXAG0AUuF9r3TVLdBrWrmQcKBZvPveC/fTTzmheWcjiJVeAKkr9PPCrU1L14q2NbLjlQgCsQCl7/3hBt/eMJ+LMffynNEVaALj++CvYv2bvvNtu+jYQptsP5rfBdPvB/DaYaP/KlStbJ06cWLbzmoXHyJWZiTTtrKrs+IfQ0hbNSOXjCVWAx3H67HAriUj2YASAtdvWpURoeGgY+w7fK89WC4IgCNkwVIg6Ml0PqwoQLPG65dAW7hiesywLX3lH5FxP20Es3/h/qdeT9jgMj2XkWyMIgmAcRva2GUJUU0ZZyJ86bu4csJBDFu6EneDNNCGavOfh+TJVEARB2AlmChEdQlQ9opzyNCHqEjmXwwZ56+s/oaHdCYWvCJRzYO2++TRXEARB6AEzhSht2Gz4qKoePSJfWgh3d0KUvs/QUaMPwevx5stUQRAEYSeYKUR0CMXwkVWUhzIDFtLZ2dCcbdsZQiTDcoIgCIXFzG0g3G25LTtBRVWQslBHM7pm4O45WGHD9o1sbakHIOQPcvBIWbwqCP3Jyxve5IF3HqO+rYHa0uGce8gZsp32EMdIIUoSstvweCzKS9M8ovbOGbh7niNKj5abuPsE/F5/lzpCcWhas4yGFxcRa6zHV1lD9YzZVBzc/YJkYeDz8oY3ufutRUTizjqcutYvufstJwO+iNHQxWghKrWcvHFlwbQ5oi4eUcccUazpy4z9hqqqQ2we/YmTog8ZlhtINK1ZRt2Td2HHnM841lhH3ZPO/lEiRuby4LuPp0QoSSQe4cF3HxchGsIYOUeUpNznrBlKD1bo7BGlzxGt21qRsd/QjoY2Qmv3pKpuNCVeP4eO+koBrBZyoeHFRSkRSmLHwjS8uKibKwQTqG/9slflwtDAbCEqcbIo9Bi+XVYJbpTd240HZ91vaLeNisNGHUTQF+hni4VciTXW96pcMIOa0uG9KheGBkYLUYW7xVBG+HanoTnL8qQCFloS2dMs+SMhJo+RYbmBhK+yplflghmce8gZlHgz87OVeEs495AzimSRMBAwWogqy50prvLS7ofmoCNgoczTkvU+0ZJ2jhh9cD9YKPSV6hmzsTp5qJYvQPWM2UWySMgHU/eaxLyjZlMTcr6TtaXDmXfUbJkfGuIYHaxQVeEIUHmGR9R1O+/kPNFhoRUsDx9HLG23iIQnRuDgJspKSvvXWKFXJAMSJGpu8DF1r0lM3v0wwKzs1UL/YawQeRIxyqoc8SjrYY4IOoRo78A/KZ9wHG/oENsbWomWtLFljOZfp8wojNFCr6g4eJoIjyAMAYwVomCsGV/pCKCTR5RFiHxpIdz71+5AfX0q8564Chsby7I4co8et1MXBEEQ+hFjhSgUbcZX5uwZFCzx4bFA7fiY4758m1fPvI9AbQ1j58xm5PRpGdkVYs0NvLPpXWyciLvxtftRFawsShsEQRAEg4UoGG3CW+oMzXk8Foe1b+ArbZ+zZo9TaPeVEYy1sN+9T3EsULZ7ZnaFNzel5ZaTaDlBEISiYmzUXCjWjK+0I8DggKatrBtxDO3+crAs2v3lvF89mdcfejEjzU9Tcz1rtujU8VFjDi2o3YIgCEImxgpRukcEsLHqYBKeTAcv4fGh/Qpv2lYQq+NNxO0EAPsNH0etLKQTBEEoKsYKUSjWjLc0lDoO+7IvVm33leEtqwJ3M73Vfjt1ToblBEEQio+Rc0QWEIo24SvtEB+vJ0Hc7rqhXXmpB8vjxVtWRUvrDj5My9Q9acxhhTBXEPLCyxve5MF3H6eu9UtqQtWcd+iZshBUGBQY6REFYi34ExHq33qro2yfWmx3yC2Jzwsn/IszB+Qtr+bl6lJiHscz2qtqD3avGFk4owVhF0hun1DnJgetb2vg7rcW8fKGN4tsmSDsOkYKUUm8HYCP75zP1qXLABi+5zA+tiDshmVXVYc47ZzDmDBxDAAt5ZW8PKxjKO+0A08osNWC0Hd62j5BEEzHyKG5JIlwmE8XLGLk9GmUh/x8ic2X2Jw1c3/OPzVzS4d/+CNEYo7u7lFSybFjjyqGyYLQJ2T7BGEwY6RHlE64ztkWoKyH7ApbmrfxSqzjC3tm6Z54PMY3XRhCyPYJwmDG2N54R2AEr+51NnW7TwB63pPo4TVLiLtDduPaIoyPWIUzVBDygGyfIAxmjB2ai3n8tPvLeS84kX1WbqSsNHsG7k8aNvLqho6ghln1zcSD2wtqqyDsKsnoOImaEwYj5gqR++swFocXntKcNKdjTVD6nkQPrn6sI69cc5i92mPEmxsKa6xgJMlw6frWL6kpHc65h5xR1I5/6l6TmLrXJCIR54eWbKEgDBbMFSJPhwe0o6GNsmDXXVrf3/ohb3/+HgAWFid92QwgQmQIxRSCZLh0MlKtrvVL7n5rEYB4IYKQZ4ydI4p5On4NVlWHMueI2qPYts2idx9LlU3b6yhGReIAxFt2YCfihTNW6DXp62ZsOoSgUOtmJFxaEAqHsUIUd4XI7/cyc5bqsjneW5v+j3X1/wTA5/HxrQmn4ylNbvdgE2/ZUWiThV5QbCGQcGlBKBxGD81VVYeYOUulFq2W+DxEYgli8TgPvPtEqu5J+01nRFkN4fJhRFobAYg3fYmvQkJf+0p/D5sVWwhqSoenshh0LhcEIb8Y6xGdcu5R/ODqr6VECKDcjZzz1m5mc9MXAIR8Qf7lKyc75Wk7tcZknqjPvPrZin4fNiv2uhkJlxaEwmGsR5SeeTtJWcjPl02t+PZYnyo77cATqAyUO9eUZ26QN1jpb2/lf99b0u2wWb6ec+4hZ2QEC0BhhSA9XHqgRM0JwmDFWCG6492HmbGHN6NjKA+V4NvtUzwBJxddVaCCrx8wM3X+83gbFe7rx1f9lXE1w7p0LPnoxHflHrv6/EJEe9W3ZRfxfA6bDQQhOKypnb0+qSPWWI+vMkH13u0Fe3Z3NK1ZxpcvLiLeWIevspbqGbOpOHhasc0ShF1ip0KklPIAdwCHAmHgIq31+rTz3wXmATHgOq31EqVULfAAEAI2AxdqrVt7U3dndm2LN3PXmwtpDrewf83efNm2nc+9/4dv9MepOqP84wj6g4DTQb+99X2+nmx4WzN3Ll8AdHR6L294kzuXLyBmxwCnE+9cZ2fMf/Ep/vHFk+CJp+6RqxDk4/ndTfLfu+KveevEa0LVWcUo38Nm8brdaX9nOm0NbbRXh4iP3h32yusjuqVpzTK2/O0OPAlnKUCssY4tf7sDoGgd/0C0qbe8tPIzVj71BNMTy6n2tBAPVbP7Sd82xv7eYFpb+7uv7+nZlm3bPZ1HKfUN4HSt9QVKqSnAVVrrM9xzo4DngCOBIPCK+/q3wCqt9X1KqSvdRj2Ya12t9S092fTW35bYf/7icbZXdq+jiUiA6D+/wojAKMbWjuSdz9fyFXsjF3r+D4BGO8BnVGJhURl0hu4a25tTi18z3iQsKoMVXco7E4nGaYu1gpXlHpaHykA5yffbsrqmGWpsb9ql5wPsaG/s9lzIW0aJv+ueTb3Btm3aoxEiiXCXcyWeAKGSwC7dP0kkGqetPZbxblhAKOjLSxsg+2eQxImqzPbdsNyNFgvPQLSpN0SicaLhMAErQuY7b+EJhLB85izQ3dm/oYHY1rIZ57ROnDgx+w6i9H9fr7Xu2mm45DI0dyzwNIDW+g2l1JFp5yYBr7oPCCul1gOHuNfc4NZ5yn39US/q9ihEABF/z/niohvGk9ixG1uw2bJ1CzCc7b4YuBHclVaYg9jmHIS3Ov/v6ZbJOjujp/CP8Jaer83D88f0dI/ENuefST7I9hyb/N0fINt3NZHnZ3RHT9+McF0BDMjCQLSpt/i7KY+7f4OJAdbWHGbF+7uvf4tuyEWIKoH0RTdxpZRPax3Lcq4JqOpUnq1sZ3W7oJSaC8wFWLRoEfNGfQcAC8v2e3xtkUS0NOOCfTpelnj8rZFErBR2YwsTOt3ZosTjawVw6mT/xZms0xORWKIUT4Js90ja2eP1u/h8gISd8MXseNcuPOFx7uPz5HSfnojEEqXdncvH/Qv1jJ6wY9FSK8tnYWNh+fz9/vxsDESbekMkligtIdZ9BV/JgG9DrgzEtkYikZBSakVa0Xyt9fy04/7u67slFyFqBNLHhTyuYdnOVQDb08rbspTlUrcL7hs2P9s5QRAEYedorXs63d99fbfkso7oVeAUAHfccHXauTeBqUqpoFKqChgPrEm/BpgFvNzLuoIgCEJh6e++vltyCVZIRlIcgjMzcKH74PVa6yfc6Ii5OKJ2g9b6r0qp3YC/4ChhHXCe1rqlN3V7NEoQBEHIK/3d1/f07J0KkSAIgiD0J8am+BEEQRAGByJEgiAIQlExJsXPzlb9FsGet+kIUfwncDfwR5yVxM9qrX/Znc3uRGCf6+6i3ZOBm7TWxyml9gPuw4kZXwNcorVOKKV+AZzqPvM/tdZv9lfdPLThCOBvwDr39J1a64cHYhuUUn7gHmAcEACuA97vD5v66zPopg0bMecz8AJ/BhTOap4LceZD8m5Tf38PBhMmeURnAkGt9dHAlcDNxTJEKRUE0Fof5/5dCNwFnIezwGuy20F2Z/Ou1u2r3VcA/42z2hng98DVWuupOF/GM9z7TwcmA/8K3N7PdXe1DUcAv0/7LB4ewG34N6Devecs4LZ+tKm/PoNsbTDpMzgNQGv9VeBa9xmmfQaDDpOEKGPVL076iGJxKFCqlHpWKfWCUmoaENBaf6S1toFngK+RxWalVGUe6vaVj4BvpB1PBJa6r58CjnfteFZrbWutPwV8SqkR/Vg3H204VSm1TCn1P0qpigHchkeAa9KOY/1oU399Bt21wYjPQGv9GO7CeJzMhVv60ab+/B4MKkwSoqyrfotkSyvwO+Ak4GLgXrcsSXcrjONuWeMu1u0TbghlNK3IcgWuJzuS5f1Vd1fb8CbwE631NOBj4BcDtQ1a62atdZPbUT8KXN2PNvXLZ9BNG4z5DNw2xJRSfwH+5LbBqM9gMGKSEPW06rfQfAgsdH/VfIjzDys99XR3K4w9Wcr6UjdfJHKwI1neX3V3lcVa65XJ18DhebKrX9qglNoTeBFYoLV+oB9t6rfPIEsbjPoMALTW5wMH4MwXpW9uZsRnMNgwSYh6WvVbaL6DO4ejlBoNlAItSql9lVIWjqf0Mlls1lo3ApFdrJsv3lZKHee+Tq6KfhU4SSnlUUqNxRH8un6su6s8o5RK7m/xNWDlQG2Du/jvWeCnWut73GKjPoNu2mDSZzBHKXWVe9iKIxYrTPoMBiPGRM3h/NI6QSn1Gh2rfovF/wD3KaVewYmI+Q7OP+hFgBdnDHi5Uuotstt88a7UzWM7fgT8WSlVAqwFHtVax5VSLwOv4/xQuaSf6+4q3wNuU0pFgC+AuVrrxgHahp8B1cA1SqnkPMsPgFsN+gyyteFy4A+GfAb/H3CvUmoZTn7s/3Tvbfr3wGgks4IgCIJQVEwaUky+pAAAAEFJREFUmhMEQRAGISJEgiAIQlERIRIEQRCKigiRIAiCUFREiARBEISiIkIkCIIgFBURIkEQBKGoiBAJgiAIReX/B6odZgvzsY2DAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "SO4_t = np.array([0.013, 1.013, 2.016, 3.005, 4.985, 6.851, 8.865, 11.882, 14.899, 18.713, 23.057, 28.923, 35.955])*24*60*60\n",
    "SO4 = np.array([1.643, 1.646, 1.565, 1.678, 1.484, 1.360, 1.121, 0.754, 0.336, 0.011, 0.014, 0.003, 0.002])*1e-3\n",
    "NO3_t  = np.array([0.036, 1.054, 1.988, 3.053, 5.008, 6.859, 9.018, 11.864, 14.857, 18.734, 23.065, 28.895, 35.898])*24*60*60\n",
    "NO3 = np.array([1.710, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0])*1e-3\n",
    "Fe2_t = np.array([0.017, 0.988, 2.003, 3.033, 4.987, 6.916, 8.909, 11.923, 14.896, 18.653, 23.070, 28.984, 35.975])*24*60*60\n",
    "Fe2 = np.array([1.588, 2.307, 5.167, 10.080, 15.267, 16.464, 18.682, 18.846, 16.893, 18.819, 19.503, 19.049, 19.254])*1e-3\n",
    "CO2_t = np.array([0.156, 1.032, 2.070, 3.072, 5.069, 6.893, 8.919, 11.928, 14.948, 18.766, 23.078, 28.996, 35.988])*24*60*60\n",
    "CO2 = np.array([1.076, 2.109, 3.247, 5.008, 7.323, 8.342, 9.048, 9.648, 10.971, 10.958, 11.512, 10.448, 10.802])*1e-3\n",
    "CH4_t = np.array([0.049, 1.016, 2.018, 3.011, 5.014, 6.872, 8.894, 11.903, 14.888, 18.724, 23.030, 28.954, 35.978])*24*60*60\n",
    "CH4 = np.array([0.005, 0.002, 0.005, 0.005, 0.005, 0.005, 0.022, 0.044, 0.077, 0.091, 0.168, 0.347, 0.637])*1e-3\n",
    "\n",
    "fig, ax1 = plt.subplots()\n",
    "ax2 = ax1.twinx()\n",
    "ax1.plot(bl.time, bl.NO3.concentration[0], c=sns.color_palette(\"deep\", 10)[0], lw=3, label='NO3')\n",
    "ax1.scatter(NO3_t, NO3, c=sns.color_palette(\"deep\", 10)[0], lw=1)\n",
    "ax1.plot(bl.time, bl.SO4.concentration[0], c=sns.color_palette(\"deep\", 10)[1], lw=3, label='SO4')\n",
    "ax1.scatter(SO4_t, SO4, c=sns.color_palette(\"deep\", 10)[1], lw=1,)\n",
    "ax1.plot(bl.time, bl.CH4.concentration[0], c=sns.color_palette(\"deep\", 10)[2], lw=3, label='CH4')\n",
    "ax1.scatter(CH4_t, CH4, c=sns.color_palette(\"deep\", 10)[2], lw=1,)\n",
    "ax2.plot(bl.time, bl.Fe2.concentration[0], c=sns.color_palette(\"deep\", 10)[3], lw=3, label='Fe2')\n",
    "ax2.scatter(Fe2_t, Fe2, c=sns.color_palette(\"deep\", 10)[3], lw=1,)\n",
    "ax2.plot(bl.time, bl.CO2.concentration[0], c=sns.color_palette(\"deep\", 10)[4], lw=3, label='CO2')\n",
    "ax2.scatter(CO2_t, CO2, c=sns.color_palette(\"deep\", 10)[4], lw=1,)\n",
    "ax2.grid(False)\n",
    "ax1.grid(lw=0.2)\n",
    "ax1.set_ylim(0, 2.5e-3)\n",
    "ax1.set_xlim(0, 40*24*60*60)\n",
    "ax2.set_ylim(0, 25e-3)\n",
    "ax1.legend(frameon=1, loc=2)\n",
    "ax2.legend(frameon=1, loc=1)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-07-02T18:08:38.356342Z",
     "start_time": "2020-07-02T18:08:38.278358Z"
    }
   },
   "outputs": [],
   "source": [
    "from porousmedialab.calibrator import Calibrator\n",
    "from porousmedialab.metrics import rmse, norm_rmse\n",
    "# available metrics:\n",
    "# percentage_deviation, pc_bias, apb, rmse, norm_rmse, mae, bias, NS, likelihood, correlation, index_agreement, squared_error, coefficient_of_determination, rsquared"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-07-02T18:08:38.359643Z",
     "start_time": "2020-07-02T18:08:38.357575Z"
    }
   },
   "outputs": [],
   "source": [
    "calibrator = Calibrator(bl)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-07-02T18:08:38.364685Z",
     "start_time": "2020-07-02T18:08:38.361060Z"
    }
   },
   "outputs": [],
   "source": [
    "calibrator.add_measurement(name='SO4', values=SO4, time=SO4_t)\n",
    "calibrator.add_measurement(name='NO3', values=NO3, time=NO3_t)\n",
    "calibrator.add_measurement(name='Fe2', values=Fe2, time=Fe2_t)\n",
    "calibrator.add_measurement(name='CH4', values=CH4, time=CH4_t)\n",
    "calibrator.add_measurement(name='CO2', values=CO2, time=CO2_t)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-07-02T18:08:38.368238Z",
     "start_time": "2020-07-02T18:08:38.365930Z"
    }
   },
   "outputs": [],
   "source": [
    "calibrator.add_parameter(name='k1', lower_boundary=1e-8, upper_boundary=1e-3)\n",
    "calibrator.add_parameter(name='SA', lower_boundary=1e-8, upper_boundary=1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-07-02T18:08:38.373741Z",
     "start_time": "2020-07-02T18:08:38.369611Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "::::: norm_rmse =  6.2062e+00\n"
     ]
    }
   ],
   "source": [
    "calibrator.estimate_error(metric_fun=norm_rmse)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-07-02T18:08:58.850673Z",
     "start_time": "2020-07-02T18:08:38.375238Z"
    },
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Converged (|f_n-f_(n-1)| ~= 0)\n",
      "::::: norm_rmse =  1.4307e+00\n",
      "Calibrated parameters:\n",
      "\tk1 =  1.1160e-06\n",
      "\tSA =  6.0000e+02\n"
     ]
    }
   ],
   "source": [
    "calibrator.run(verbose=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-07-02T18:08:59.221196Z",
     "start_time": "2020-07-02T18:08:58.852566Z"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "'c' argument looks like a single numeric RGB or RGBA sequence, which should be avoided as value-mapping will have precedence in case its length matches with 'x' & 'y'.  Please use a 2-D array with a single row if you really want to specify the same RGB or RGBA value for all points.\n",
      "'c' argument looks like a single numeric RGB or RGBA sequence, which should be avoided as value-mapping will have precedence in case its length matches with 'x' & 'y'.  Please use a 2-D array with a single row if you really want to specify the same RGB or RGBA value for all points.\n",
      "'c' argument looks like a single numeric RGB or RGBA sequence, which should be avoided as value-mapping will have precedence in case its length matches with 'x' & 'y'.  Please use a 2-D array with a single row if you really want to specify the same RGB or RGBA value for all points.\n",
      "'c' argument looks like a single numeric RGB or RGBA sequence, which should be avoided as value-mapping will have precedence in case its length matches with 'x' & 'y'.  Please use a 2-D array with a single row if you really want to specify the same RGB or RGBA value for all points.\n",
      "'c' argument looks like a single numeric RGB or RGBA sequence, which should be avoided as value-mapping will have precedence in case its length matches with 'x' & 'y'.  Please use a 2-D array with a single row if you really want to specify the same RGB or RGBA value for all points.\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7ff748c99510>"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaIAAAD7CAYAAAAo/ZDkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOydeXzUxf3/n3sk2dzXJlzhRgaUQ0A5lLuIoqK2ihelYrVg61dttVpt1f7aetS2autVi61VAcGqRRAVxQvkVqDcDJcCIZBkc23OPT+/P3az2Zxskk32yDwfrvl85jOfmfcnS+b1mZn3vEenaRoKhUKhUIQKfagNUCgUCkXXRgmRQqFQKEKKEiKFQqFQhBQlRAqFQqEIKUqIFAqFQhFSlBApFAqFIqQYz5ZBCKEHXgJGAjbgdinlEb/rPwEWAk7gMSnlaiGEGXgTiAfygFullFVCiF8AN3pv/VBK+TshhA7IBQ570zdLKR8KzuMpFAqFIhDa2Nb3AV7FoyU6YIGUUgoh7gVuAwq9ty+UUsrm6j6rEAHXACYp5QQhxHjgaeBqr2HdgbuBCwATsEEIsRZ4FHhTSvmaEOJBYKEQYiUwFxgHaMBXQogVQBWwQ0o5OwBbFAqFQtExtKWt/wPwgpTyPSHEpcCTwA+A0cCPpJTbA6k4kKG5icAaACnlFq8htYwFNkopbVLKMuAIMML/HuAjYAZwErhMSumSUrqBGKAGGAP0EkJ8IYT4UAghAjFcoVAoFEGlLW39fcAH3jxGPG06eNr1h4QQG4QQZx3hCqRHlAKU+Z27hBBGKaWziWvlQGqD9HIgVUrpACzeobg/AzullIe8SvuklPJtIcREYAlwYUMjhBALgAUAS5cuHaPXq+kthUKhCBSbzabNnz9/h1/SIinlIr/zVrf1UkoLgLcD8Rc8vSqA5cCLgBVYIYS4Ukq5ujnbAhEiK5Dsd673GtbUtWSg1C+92i8NIYQJz3hiOfAz7z3f4BlzREq5QQjRSwihk1LWiz3k/YUtAtixY4c2evToAEwPP+x2OwCxsbEhtqRtRLr9EPnPEOn2Q+Q/QyTav3379mop5QUtZGlLW48QYhqeuaV53vkhHfBXb88JIcQHwCigWSEKpFuxEbjcW+B4YI/ftW3AJCGESQiRCgwF9vrfA8zCMx+kA1YCu6SUC6WULu/13wI/95Y/EjjRUIQUCoVC0eG0uq33itDf8Ey7fOPNm+K9luRt96cDLc4V6c4W9NTPk2IEHq+IW73GHpFSrvJ6UizAI2pPSCnfFUJ0A17Ho5oW4GZgJrAM2OJX/EPAQTzDcUl4ekZ3SikPtmST6hGFjki3HyL/GSLdfoj8Z4hE+7dv3141ZsyYxOaut7Gt3wXEAWe8xUgp5UIhxDw8zg024DMp5W9bsu2sQhSOKCEKHZFuP0T+M0S6/RD5zxCJ9p9NiEJJIHNEEYHD4SA3N5eampqzZw4htcKv0+kCvsdkMpGTk0NMTExHmaVQdCna21605e+4s4jE9iJqhCg3N5fk5GT69esXlv84anG73QAE6vWnaRpFRUXk5ubSv3//jjRNoegytLe9aO3fcWcRqe1FeP0W20FNTQ2ZmZlhLUJtQafTkZmZGfY9PYUiklDtRXgRNUIE4dlNDgbR+lwKRSiJ1r+rSHyuqBIihUKhUEQeUTNHFGq2bt3KnXfeyfvvv0+PHj0A+Mtf/sKAAQO49NJLefbZZzlw4AA6nY7ExEQefPBB+vfvj8vl4uGHH+bbb7/FYDDw5JNP0qdPnxA/jUKh6Ghyc3O56qqrOO+883xp48aN4//+7/8a5T1w4AB/+MMfMBgMxMbG8tRTT2E2mzvT3A5FCVEQiYmJ4aGHHuLf//53ve7xI488wqhRo3j44Ydxu90cPHiQO++8k7feeoutW7cCsHz5crZu3cqTTz7J3//+91A9gkKh6EQGDRrE4sWLz5rv8ccf55FHHmHo0KEsX76cV155hYceip5NCqJSiFZ8eYRlnxyk2uY6e+YAiY8zcNPMIXx/6qBm84wfPx63283SpUv54Q9/CEBJSQmHDh3imWee8eUbMmQI06ZN45NPPuHaa69l6tSpAOTl5UXVW45CEQmcem8VJ5a9hTuIE/x6k4k+N91Ar2uuavW9Tz/9NF9//TWapjF//nxmzZrFM888Q3Z2NgAul4u4uLig2RoORKUQvbfuSFBFCKDa5uK9dUdaFCKA//f//h9z5sxh4sSJgMfNs3fv3o3y9e7dm7y8PACMRiO/+tWvWLt2Lc8991xQ7VYoFC1z6r1VQRUhAHdNDafeW3VWITpy5Ajz5s3znc+ZM4fc3FyWL1+OzWbj+uuv5+KLL/aJ0I4dO1iyZAlLly4Nqr2hJiqdFa6ZMoj4OENQy4yPM3DNlJZFCCA9PZ1f//rXPPjgg7jdbhwOh09w/Dl+/LhvLgngqaee4uOPP+aRRx6hqqoqqLYrFIrm6XXNVehNpqCWqTeZAuoN1Q7N1X7y8/PZt28f8+bN4/bbb8fpdPrajw8//JDf/va3LFq0iIyMjKDaG2qiskf0/amDztpz6UimT5/O2rVrWbFiBffffz99+vRh6dKlzJ07F4B9+/bx+eef89Of/pT33nuP/Px8Fi5cSHx8PDqdDoMhuCKqUCiap9c1V7V6CK2jFrQOGDCAcePG8Yc//AG3281LL71ETk4OK1eu5K233mLx4sWkpaUFtc5wICqFKBz4zW9+w5YtnviuTz31FH/605+YM2cOer2elJQUXnrpJVJSUpg5cyYPPfQQc+fOxel08utf/zrqxn8VCkVgTJ8+nW3btnHzzTdTVVXFjBkziI+P5/HHH6dHjx7cddddAFx44YXcfffdIbY2eERN0NMDBw4wdOjQEFkUOG19kwqX54vEYI8NifRniHT7IfTP0N6/p3AN8VNLU88XzkFPw/O3qFAoFIougxIihUKhUIQUJUQKhUKhCClKiBQKhUIRUpQQKRQKhSKkKCFSKBQKRUhR64iCyKJFi9i0aRN6vR6dTscvfvELhg0bxkcffcSSJUvQ6/U4nU7mzJnDD37wg3r3PvLII6SmpvLLX/4yRNYrFIrO5PDhw/z5z3+murqaqqoqpkyZwl133UVJSQlPPfUUeXl5uFwuevTowYMPPkhWVhbl5eXcf//9VFRU4HA4ePDBBxk1alSoH6XdKCEKEkeOHOHzzz9n2bJl6HQ6Dhw4wK9+9SseeOABli9fzssvv0xycjJVVVXcc889xMfHM2vWLMATefvQoUNceOGFIX4KhULRGVitVu69916ef/55+vXrh8vl4p577mHZsmWsXr2aH//4x8yYMQOATZs2sXDhQt5++23+/e9/M378eObPn8+xY8e47777WLFiRYifpv1EpRCVbllFyVdvodmDF8hQF2sifdINpI1vOhRIRkYGeXl5vPPOO0yePJmhQ4fyzjvvcNddd/HLX/6S5ORkAEwmEw888AC/+93vmDVrFjt37mTXrl3ccMMNHDt2LGj2KhSKwNj85VHWfXIIexADJcfGGZgyczATpg5s8vpnn33GuHHj6NevHwAGg4GnnnqKo0ePsm7dOp8IAVx00UX06dOHr7/+mvnz5/sWAUdTFO6onCMq27oqqCIEoNlrKNu6qtnrGRkZ/P3vf2fHjh3ccMMNXHbZZXzxxRecPHmy0UZ3tZG3CwoKeOGFF3j00UeDaqtCoQiczeuOBVWEAOw2F5vXNf9iWVBQ0Cgqf2JiIrm5uS1G609JScFkMlFYWMj999/PvffeG1S7Q0VU9ohSx13VIT2i1HHNB0Y8fvw4SUlJPPnkkwDs2bOHBQsWMGTIEE6dOkVqaqov73fffUePHj1Ys2YNJSUlLFiwgMLCQmpqahgwYECj+SOFQtFxTJgyoEN6RBOmDGj2es+ePdm/f3+9tJMnT2I2mzl16lSj/MePH+eiiy4CQErJvffeywMPPMDYsWODZnNI0TQt4j7bt2/XGrJ///5GaZ3Jxx9/rM2fP1+rqanRNE3TysvLtUsuuURbu3atdsstt2jl5eWapmma1WrVbrvtNu2DDz6od/+7776r/fnPf262/FA/Xy02m02z2WyhNqNdRPozRLr9mhb6Z2jv35PL5dJcLleb7y8vL9euuOIK7fjx45qmaZrdbtd+9rOfaUuWLNGuu+467bPPPvPlXbdunXbVVVdpTqdTO3z4sHbppZdqBw4caLH8pp7vm2++qdTCoP1u6hOVPaJQMHPmTI4ePcqcOXNISEhA0zQeeOABZsyYQVVVFbfffjs6nQ6Xy8V1113H5ZdfHmqTFQpFiEhKSuKPf/wjDz/8MJqmUVlZybRp07j55pu57LLLeOKJJ/jHP/4BQPfu3Vm0aBEGg4Gnn34au93O448/7ivn73//eygfJSio6NudjIq+HXoi/Rki3X4I/TOo6NvhRXj+FhUKhULRZVBCpFAoFIqQooRIoVB0SSJxWiIQIvG5lBApFIouh8lkoqioKCIb7ZbQNI2ioiJMJlOoTWkVymtOoVB0OXJycsjNzaWwsLBN99cKmE6nC6ZZQcFkMpGTkxNqM1qFEiKFQtHliImJoX///m2+P9Ref9GGGpoLIocPH2bBggXMmzePa6+9lueee46TJ09y/fXX18u3fPlynn/+ed+52+3m9ttvZ9myZZ1tskKhUIQc1SMKEs1F092wYcNZ7/3rX/9KWVlZJ1ipUCgU4UdUCtH7Bz/l7X2rqXHaglamyRjHnPOuZPaQGU1eby6abkFBQYth2tesWYNOp2Py5MlBs1WhUCgiibMKkRBCD7wEjARswO1SyiN+138CLAScwGNSytVCCDPwJhAP5AG3SimrhBC/AG703vqhlPJ3Qoh4YAmQDZQDt0gp2zaD6GW1/DSoIgRQ47SxWn7arBA1F003JiaGI0eOMG/ePMAzyVlQUMDs2bM5dOgQq1ev5rnnnuPFF18Mqr0KhULRGtrY1vcBXsWjJTpggZRSCiFmA496874qpXylpboDmSO6BjBJKScADwJP+xnWHbgbuBi4FHhSCBHnNeBNKeUkYCewUAgxAJgLXARMAGYKIUYAPwX2ePO+ATwcgE0tcqWYgckY3H06TMY4rhRNixB4oumeOXOmXtrJkyc5ffo0gwYNYvHixSxevJg33niD+fPnA/Dee++Rn5/PLbfcwooVK3jttddYv359UO1WKBSKAGlLW/8H4AUp5VTgCW96DPAsMBOYAizw3t8sgQzNTQTWAEgptwghLvC7NhbYKKW0ATYhxBFghPeeJ7x5PvIevwBcJqV0eR8sBqjx5v2TX95HArDJ57VSi6ZpvvhPVwyezhWDpwdSTKupraMhU6ZM4eWXX+aGG26gT58+OBwOnnzySS666KJ6tvlHnPXfFvyFF17AbDYzceLEJuvQNK3RM4cCh8MRahPaTaQ/Q6TbD5H/DJFufzO0pa2/D6id4DbiadOHAkeklCUAQogNwCTg7eYqDkSIUvwqAnAJIYxSSmcT18qB1Abp5UCqlNIBWIQQOuDPwE4p5SEhRKO8TRkhhFgALABYsmRJAGZ3LrXRdB999FHcbrcvmu7kyZNZuXJlqM1TKBRdHKvVahRCfOOXtEhKucjvvNVtvZTSAiCEEMBf8PSqsprK25JtgQiRFUj2O9d7DWvqWjJQ6pde7ZeGEMKEZzyxHPhZE2X48jbE+wtbBJ7o2w3993U6Xcgj4Q4fPpw33nijUfrbb9e9CLjdbm666aZGtt59990tlq3T6cJqzUI42dJWIvEZCtat5/gbS7EXWYgzm+kzby7ZUyLX0SUSvwN/Isn+lJQUp5TyghaytKWtRwgxDc/c0jzv/FBcc3mbI5CWeyNwubfC8cAev2vbgElCCJMQIhVPl2yv/z3ALOArb09oJbBLSrmwdoiuqbwB2KRQdDkK1q3n6IsvY7dYQANboYWjL75MwTo1r6gICq1u670i9Dc80y61va0DwDlCiAwhRCwwGdjcUsWB9IhWAJcIITbh8Yq4VQhxL54xwFVCiOfwiIce+I2UskYI8RjwutfLwgLcjKfLNgWIE0LM8pb9EPB3b94NgN2bV6FQNODE4qW4bfW9Qd02GycWL43oXlE0U7BuPScWL8VmKSLOnBnuPdi2tPV/BWLxtOEAUkq50Hvfx968r0opG+9/7kdUbYw3ZMiQsIz95E9bNtTSNI2DBw+qjfGCRKQ+w8ZrroOm/l51Oi5+753ON6gdROp3UEsg9tf2YP1fHvRxcQy8846QiJHaGK8TUNF0FdFOnDmzVenhSsG69ez66V18ff1NfHP7wqgdWmypB6uoT9REVmhvNN3Ooi1ReyMxmq4i+PSZN7fJN+w+8+aG0KrW0bCXUDvPBYTzkFWbsFmKWpXelYkaIWpvNN3OItKHJGqJsLHvNhFuz1hbt8drLjxsai1daZ4rzpyJrdDSZLqiPlEjRIo6OroBtXy1geP/+GdUv9WG65t79pTJpE0YD0Tmy0xX6iVEQw+2s4iaOSKFh9oG1FZoAU3rEBffU2++FfVj32p8v2OIlnmuQMieMpmBd95BXJYZdDrisswhc1QId1SPKMrojKEPe1Hj4QaIrrfarvTm3pl0tV5C9pTJSngCQPWIoozOaEBjM81NpkfTW21XenPvTGp7CbFm1UtQ1KF6RFFGZ0yQ9rr5hnpzRBB9b7Vd7c29M4n0eS5F8FE9oiijz7y56OPqb4ER7AbUPGli1I99q/F9haLzUD2iKKO2oexot+OuMPbdFZ5RoQgHlBBFIaoBVSgUkYQamlMoFApFSFE9IkXYo7nduB0ONIcTt9OJ5nSiOR2+Y7fDm+Z2o7lc4HZ7jv3PXS40twaahsNhB03DaDB409ye0Etuz8651H7QvD8855o3zWMUvgCkvviG/nEOG8Q8bDJPkw/bzH1+uFyeLWIMBmPTZQYh3mJHx2x0u9xoaBgMhg6tp6NwuTy72ARkf7jEvxwqQm1BsyghUnQ4mtuNs6ISR1kZDmsZjjIrjjIrTqvnp8Nahqu6BldNDW6bre5nteenOwy2SFcoIh3To78OtQnNooRIETSclZVU556i6mQu1bm5VOWeojo3l5r8AvBuf6FQKBQNUUKkaBNupxPr/gOUbN9B5bFvqTqZi6OkpMPq08fGoosxojca0Xk/dccx6I0GdAYD6PXo9Hp0BgM6vd5z7j32nOtwuzXP1vJGoycKul6HTqf3bAWm06PTgfd/nv9817xpvkOd9wBvWv1zvwt115uiYfpZIrO7XW7QNRgWamUZoaZ2aMtojMwmqFVDcxAW30dBqA1ogcj8V6AICQ5rOZZt2yj9ZgfWXbtxVVUFfK8hMYGYlBRiUlKJSU3BmJJCTGoKMampxKQkY0hIxBBvwhAXh95kwmCKQx/n/Rkb6xGRIBHpEdAj3X6I/GeIRPsLtm8PtQnNEpFCZHO4+PFjn/CjWUOZOqZ3qM2JaqpOnKT4628o+WY71oOyxSE2XUwM8T17EJ/Ti4ScHOJzckjonYOpZw8MDRbZKhQKRS0RKUQAhSXVvPD2LgAlRh1AuTzE8SVvUrZ7T7N54rKzyLjwAlJHjiChT29M2dme4TGFQqFoBREpRAbcfD/hawBOf/QNluK+9a7r4xIwppiJSc3CmJqFMcWMzhgTClMjjsrvjnNi6TKKt33d+KJOR9Lgc0gdM4qs8eNJ6NO7VTvNtpZw25hOoVB0DBErRFNNB3zn1m3Nv7V70GFISvOJUkxGT5JHTiMmvXujnOV711PyxVKc1iKMKZmkT5tL8rDob/yqT5/mxJtvYflqQ/11D3o9mePGkjH2QtLHjEKLjwc6fmw8XDemUygUwScihaj1aLgqSnBVlGA7dQiA0s0rSBl1CWkT52BMSgM8ImT54GU0p6fxc1otWD7wNH7RKkY2SxEn//M2+Ws/qz//o9NhnjSRPjddT3zPnr5keyet6elKW0orFF2diBQiF3pWVF6A0aBn4vk9GZiT5ndVw11VjrOsEEdZIc6yQlwVJaA1mGR3u7BuX0P57i9JHTebtPFXUfLFUp8I+Upz2ij5YmnUCZHmcnHizeWcWvk+msNR71r6hRfQ94c3kdivX2iMQ21Mp1B0JSJSiAy4mZZwEOfIaxgz+/u+9OaG1TSXE2d5Mc6yQpxlBZT/7zNqTnqG9jRHDaUb3sa642PcVdYm63NaO6/x64x5EUd5OfLPz1C2a3e99JRh59F33lxShoQ+FEhn7KukUCjCg4gUIoA0XQW6/W9T3j+T5GGTzzqsFpOWTUxaNnAeScOnUn1kB8VfLsFecAKgWRECMKZ0TuPXGfMiVSdzOfD4k9ScPuNLSzpnEH1/eDOpI0d0qPNBa1Ab0ykUXYeIFSKoP2zWmmE1nU5HwjljiB94PhV7v6Jk3TKc1sZv3wA6Yxzp0zqn8evoeZGS7TuQf3m23kLU3jfdQO/rrwvqgtFg0Fn7KikUitAT0UIEdcNmzQ2ftTSsptMbSB4xlcRzL8K6/WNKN76Lu7q8Xp7Y7v2J7zs8eAa3QEfNi2iaRt7K9/nu9cU+hwR9XBzn3HMX5osntKvsjkTtq6RQdA3C6zW4DdQOmzU3fBbIsJreGEvauNn0+dmLpF18LTpjnWuyLfcgJ1++i9ItK9FcjhZKaT/NzX+0Z17Ebbdz5LkX+O7fr/tEKNZsZvgfHw9rEVIoFF2HiBYi/2Gz9Glz0Rnjmr0eCHpTIhlTbybnjr+ROGS8L12zV1P82RvkvnIfVcd2Bcf4Jugzby76BqFw2jMvYi8pYe/Dv6Xg8y99aclDBCOffoqkAf3bY6pCoVAEjYgdmjOmmOstNq39GYzFqDGp2XS79n6qvt1F0Sev4rDkAuAoOsWZZb8nQYwjc8Z8r/ND8AjmvEhVbi77Hv099qK6Yb3s701n4E8XoI9RUSYUCkX4oOvonRg7gh07dmijR49u9X1tiZqguZyUffMhJev/g2av9qXrjLGkT7mR1LFXotMHHl+tM6L22ktL2X3/Q9gKvIHf9Xr6zf8RPa+6st1ecZEYdbghkf4MkW4/RP4zRKL927dvrxozZkxiqO1oiogemmsNte7dHu84zefeXb53fYv36QxG0sZdRe+fPk/S8Km+dM1pp/izN8h7/TfYC092rPGtwGWzceCxP/pESG8yce4jv6bX1bPDxjVboVAo/OkyQtSSe3cgGJPSyb7qLnre8jix2f186ba8w+T+65eUbHwXze0KpsmtRnO7Ofzsc1QcPuxJ0OsRv/wF6aNHhdQuhUKhaImInSNqLW1x724KU84Qev34KUo3raBkwzvgdoLLScmXb1J5cAv6rPHkrfokJGtfji9eStHmLb7zAbffSsaFF3RK3QqFQtFWuowQGVMym1y02paoCTqDkfRJc0gUYylc/SK200cBKN17jLLvjoE3rF1nRow+8/EnnPrve77zHrOvoMcVl3donQqFQhEMzipEQgg98BIwErABt0spj/hd/wmwEHACj0kpVwshzMCbQDyQB9wqpazy5s8CNgHDpZQ1QggdkAt4x5PYLKV8KFgPWEv6tLn1QgBB+6MmxGb3pef8JynbsoqS9W9RnuvwiVAtnRExumTn/zj68iu+8/QLL6D/rbd0WH0KhSL6aEtb73ft50B3KeWD3vN7gduAQm+WhVJK2VzdgfSIrgFMUsoJQojxwNPA1d7KugN3AxcAJmCDEGIt8CjwppTyNSHEg17jnxVCXAr8EejmV/5AYIeUcnYAtrSZYLp3+6PTG0i76PskDL6QvM33NJmnIyNGVx4/gfzT077FqokDByDu+7naKVWhULSWtrT1euAVYBzwrl9Zo4EfSSm3B1JxIEI0EVgDIKXcIoTwn3QYC2yUUtoAmxDiCDDCe88T3jwfeY+fxdNfmAH4GzcG6CWE+AKoBn7RknICOJzuNu2LEzd4PN0Hj6+XFrT9dVKyiTWbsVuaGP5LiKGmvBR9XAIOR/CiMzhKStn/+8d9seNiMjMY9MB9uAwGXB20b1Aw7Q8Vkf4MkW4/RP4zRLr9zdCWtv4I8AbwKTDEL/8Y4CGvgH0gpXyypYoD8ZpLAcr8zl1CCGMz18qB1AbptWlIKddKKRt2D04DT0opp+ERrCVNGSGEWCCE+EYI8Y1b0yix1gRgeudinjwKXcPfqB4Su9nJf+PX2PO/DVpdLpuNQ0/92Sd8epOJwQ89QGxGRtDqUCgU0YPVajXWtqHez4IGWVrd1kspS6SUnzRR3XLgDmA6MFEIcWVLtgXSI7ICyX7neimls5lryUCpX3q1X1pzfINnzBEp5QYhRC8hhE5KWW+lrZRyEbAIYPPWr7WyKhfdzOG1mMxg3UlKPyjPBbcd9LGQnAMJZnCV5lOw9FHSps4jcdQl7VoIp2kax/76PFVHj3kS9HqG/OqXpJ1zTnAeJAAiaSFfc0T6M0S6/RD5zxBJ9qekpDillC250balrW+Ed97/r1LKMu/5B8AoYHVT+SGwHtFG4HJvgeOBPX7XtgGThBAmIUQqMBTY638PMAv4qoXyfwv83Fv+SOBEQxFqCktp9dmydDpOaxEJZuh2PvQY6/mZYPbL4HJS+tm/KV71N9w1lW2ux/LVRoo2bvKdD1hwm1orpFAo2ktb2vqmSAH2CiGSvKI0nfrTMY0IRIhWADVCiE145nl+IYS4VwhxlZTyDPAcHqH5HPiNlLIGeAy4UQixEZgAvNBC+X8Epggh1gHPAPMDsAlLWfgJUXOu4IbEdGK71QUZrT60ldx/3e9z+24N9tIyji36p++828wZ9Jh1WeuNVSgUivq0pa1vhLcn9GvgC2/+fVLKD1uqOCJjzW3e+rUmCxOYf+V5oTalHg13iQWPi7j5ijtIHDKe4k9fx7p9Td0NBiPZs+8i6byJAddx8E9/oWjjZsCzncOo55/FmJAQtGc4G5EYY6shkf4MkW4/RP4zRKL9KtZcB2ApDT9nheRhkzFfcQfGFDOgw5hixnzFHSQPm4zeGIv5sp+QMftudLHxnhtcTgree5aSDe8QyAuBZdNmnwgBDLrzjk4VIYVCoegIIjayQjgOzYFHjFpam5QwZAKx3fpT9N7Tvu0lStYtw1FyhqzLF6IzNL1Fg8NazjG/RavZM6areSGFQhEVRGyPqChMhZ12VJkAACAASURBVCgQjOnd6fmjxzH1q9uCvGL3F5xe9hiu6oom7/n2n6/iKPN4T8ZmZND/1vmdYapCoVB0OBErRJbSmoCGs8IVQ3wSPW78DUkjpvnSao7vJe/1X+MoOVMvb/G2rylcV7ddxcCfLcSYFJZDvQqFQtFqIlaInC431sqOiR7QWegMMWRdeSfpU2/2pTmKTnHqtYeoOXUIAGdFJUde+ofvetbUySqitkKhiCoido4IPGuJUpPiQm1Gu9DpdKRffC0xadkUvP8CuJy4q6zkvfYQ+vhkqsp74ygpASAmLY3+t/04xBYrFApFcInYHhFAUVn4ec61laTzJpE24Zp6aVWnyynZsd93PvCOnxCTktzwVoVCoYhoIlqIwtVzrq1U7P7Sd+x2Qdl3ddcyL76IzAnjG92jUCgUkU5kC1EYhvlpD/67xZaf8MSrA9AbIWtc/2buUigUisgmooUomobmoC5EkL0cqgrr0lP6QtmGpZR981GILFMoFIqOI6KFKNp6ROnT5oIhFuvJurS4NDB5d3Yo+vifSowUCkXUEdFCFMmLWpsiedhk4gZeiqN2TasO+t5yM6acwb48SowUCkW0EdFCZCmL7EWtDdFcLgrW/8933uPyy8mafi09bnyYuF5KjBQKRXQS0UJks7uorI6eLXsL162n6oRnXE5vMpFz/XXe40QlRgqFImqJSCHS6XS+Y0uUOCy47XZOvLncd97rmquITUv1nTcnRvW2lVAoFIoIJDKFyO84WhwWTn/0MbZCCwAxqSn0vPqqRnmaEiPLmn9SKbd1mp0KhUIRbCJTiPyUKBocFpyVleS+/a7vPGfOdRgT4pvM6xOjnud4UzQKVv61Tbu9KhQKRTgQkULk3ycKxw3yWsup91bhLC8HIC47m+6XzWwxv96USPfrH8KY1g0AzWHjzFtP4CwrbPE+hUKhCEciUoiiqUdkLy0lb+X7vvM+c29EH9P05nj+GBJT6X7jb9CbkgBwVZZy5j9P4LZVdZitCoVC0RFEgRBFdo/o5Ftv47bZAEjo15esyZMCvjc2sxfdrrvfEwMIsBecIP+/z6C5XR1iq0KhUHQEkSlE/kNzEdwjqj59hvyP1/rO+86bi07fuq8kvu8wsq64o67MYzsp+vhfUbW+SqFQRDeRKUT+PaII9po78eYyNJen95Jy3rmkjxndpnKSR0wj7eLrfOfWHR9Ttm11UGxUKBSKjiZiN8aLMepxON1U1jipqnGQYDr7vEo4UXHsGJb1G3zn/W6ZV299VGtJn3IjjpLTVO7fCEDxp68Tk9aNRDG23bZ2BTRNw1bjpKbaQXWVg5pqR73j6moHtmoHbreG5tbQNM89no/n2O1u3Att3DFtf0/V5XIDoG9l7zmccLsj+xki0f6Bw9vevnQ0EStEmakmzhR5JuaLymoiToiOv7HUd5wxbizJYnALuc+OTqcja/b/4bRasOVKat26e877A3E9BrbT2ujBYXeRn2el2FJFWUkNRQWVFBVWUFRYia3GGWrzFIoOY+DwHqE2oVkiWIji/YSomt7dImfnUuvefZTu9MaU0+vpO29uUMrVG2Ppft2vOPXaQzhL831u3b1+/CffFhNdBU3TsJbWcDq3lLyTZZzOLaMwvxxrFLj7KxTRRsQKkTm1bsFnpK0lyvNbvJo9fSoJvXOCVrYhMZXk0TMp+XwJoOGqLOX0m78jZ8Gz6PSGoNUTblSU28g7WUreyVJOnywjL7eMynJbq8qIiTUQnxBDfHwMpoQYTPF1n/iEGEymGAxGPTqdpweq0+nQ6XV+5zQ5vNqOEdd6fPjuXqoq7Y3SExJjueK64cGppJNwOD29zxhjZDZBTq/9xgiyv8qRF2oTmiVyfosNMKeZfMeRtJao8uhRyvcfAEBnMNDnxuuDWn753vWUrv8P/nMRjqJT5L/7F7rP+VVQ6woVmqZRUlTFiWPFnDhWzPFjRZQUBbZ+SqfXkZYeT0ZWAlndUjBnJ5KRlYQ5K5HE5Lh2zdN1NO+8saPJ9KpKO0NHhO+wS1PY7R5BjY2NDbElbSMS7d++XQlR0Mn07xFF0FqiM+9/4Ds2T7yYuKysoJZf8sVSNGfjnkDVoW1UfbuLhP4jg1pfR7Nney6ffXgQa2kNpoQYzFmJlJVUU249e28nNs5Ij5xUevZOpWfvNLr1SCE9MwGX2/M2G0gjsmd7Lp9/JCkrqSY1PZ7pswTDxwSvB9saUtPjKStp/NKVmt50OCiFIlKIWCHy7xFFSuBTW2EhxZu3+s57XjM76HU4rUXNXitc+Td63fY0xuT0oNcbbCqsNaz7+BA7tp70rYmqqXKQe7y0yfxGo94rOmn08ApPpjkRnb5xD8fVeHSrSfZsz2X123twODwu9mUl1ax+ew9ASMRo+ixRzx6AmBgD02eJTrdFoQgmEStE/j2iSBmay3v/A/C6faaOGE7SgAFBr8OYkonTamnymquyjIKVf6XHzY+G3XyR0+nixLESjsoCjslC8k+Xt5g/zmSkd/8M+g7IoM+ATHrmpGIwBteV9vOPZL1GH8DhcPH5RzIkQlRbZ20PMdQ9NMXZCacedTgTsUJkTossZwVnZSX5n3zqO+95dfB7QwDp0+Zi+eDl+sNz+hhwezYQrDm+l9IN75I+ObhzU22hvKwGuS+fQ/vz+e6IBafDHfC99//hUvRN9HaCSVPDYC2ldwbDx+QghmcDkTU/0RUJtx51OBOxQpSaFIdBr8Pl1iivsmNzuIiLCa+3fH/y136Gq9rTgJl69SJ99KgOqSd52GTAM1fktBZhTMkkfdpcHEV5lG5423Ptq/9g6jOU+H6d62mlaRqW/AoO7j2D3JdP3ommh9kA9HqPR5rL2VicUtPjO1yEautRczKKthJuPepwJmKFyKDXkZFqotDbUBSVVdPTnBRiq5rG7XSSt6ou5E732Ze3OqZca0geNtknSLVobhc1J/dTc3wfoFHw3l/pdfvTGJPSOswOALdbI/e7EuS+M8i9Zyi2NO/dlmFOZKDIYoDIot/ATA7tOxPSORE1J9Mx+DugRPNwVTj2qMOViBUi8Kwl8glRaU3YClHRxs3YizxOBMbUFDInTex0G3R6A9lX/5zcf96Hu8qKq7KUwlV/o/uNDwd9vkjTNE6dKGXfzjz27cqjohkPN51eR98BGQwZ1p1zzs0mPTOx3vXaxilUY+yhrj8a6UrDVapHHTgRLUSZqX6ec2HqsKBpGqdWrvKdZ186E32IxvaNyRlkX30PZ5Y9BmhUf7ub0k0rSJ943VnvPRuaplFwupy9O0+x7395lBY3/X3ExhkYKLIRw7pxztBs4hNa/l0MH5MT0gYq1PVHG11puEr1qAMnooWovsNCcISoYN16Tixeis1SRJw5kz7z5pI9ZfLZb2wG6959VB49BoA+NpbsmZcExc62kjDgfNIu+j6lm/4LQMn6tzD1Hkp83/PaVF6xpdIjPjvzKMyvaDJPYlIsYlh3xLBu9B9kxhjGc3mKjqUrDVepHnXgnFWIhBB64CVgJGADbpdSHvG7/hNgIeAEHpNSrhZCmIE3gXggD7hVSlnlzZ8FbAKGSylrhBDxwBIgGygHbpFSBrTndX0X7vZ7zhWsW8/RF1/2bVRnK7Rw9MWXAdosRvV6Q9OnEpOa0m4720v6lBupOXmAmpMHQHNTuOo5chY8iz4uIaD77XYnB3bls3PbSU4cK24yT5zJyNARPRg2qif9BmaiN0ROlGJFx9HVhqsiqUfdlrbe79rPge5Syge957OBR715X5VSvtJS3YG0DtcAJinlBOBB4Gm/yrsDdwMXA5cCTwoh4rwGvCmlnATs9BqPEOJS4BOgm1/5PwX2ePO+ATwcgE1A8Be1nli81CdCtbhtNk4sXtrMHS1TdTKXkq+3e050Onpe1TEu261FpzeQfc0v0Md75tScVgtFn77e4j2appF7vISP/ruPF55Yx8rluxqJUEysgfPO78kNt17Afb+7hKtuGMmAwVlKhBQ+ps8SxDToEavhqrCh1W29ECJeCLEEuNMvbwzwLDATmAIs8N7fLIEMzU0E1gBIKbcIIS7wuzYW2CiltAE2IcQRYIT3nie8eT7yHj8LuIEZwPYG5f/JL+8jAdiE3W4nNaHOfEtptS/+U1uxWZpeCGqzFLWp7NwVK33HaWNGY8gy43A42mxfUDElk/a9Wyle/TwA5f/7lLhBYzD1P79etqoKO3t35rF7+yks+ZWNitHpdQwUZs4d2Z1BQ7OIjfV8J263C7s9PLcsD5vvoI1Esv1ieDZO51C+/Pgw5WU2UtJMTJk5CDE8u91/v51JJH8HLdCWtv4Ing7Ep8AQb96hwBEpZQmAEGIDMAl4u7mKAxGiFKDM79wlhDBKKZ1NXCsHUhuk16YhpVzrNay58n15GyKEWAAsAFiyZAlQ31mhyNr+obnYTDP2JsQoNrP1Wyg4ysqwrP/Kd9599hXtsq0jiB8ygfhDW6k+tA2AkjWv0O3WP6GLSyDvRBnbN5/g4N583K7Gm7klJsfidLix1TgpOF3O0OHdfCKkULTEeaN6MniYJ8ZiTExk7SMWyVitVqMQ4hu/pEVSykV+561u671i84kQYn4L5TTbrtcSSMthBfw3+9F7DWvqWjJQ6pde7ZcWSPnN5vX+whYB7NixQ4uNjSU704hO59kFs6zChk5vJKYdYV76/mhuvTkiAH1cHH1/NLfVq9jPfPo5mvetKWnQQDJGjqgX2TlcVsVnX34HJ3MP4q6yYi8v5eu33kNWDiLvZFmjvDGxBoYM70ZKahxb1x/H6V1sai2tYc2KAxiNxogZD4fw+Q7aSqTbD5H/DJFkf0pKilNKeUELWdrS1gdSztk0ICAh2gjMBv4jhBgP7PG7tg14XAhhAuLwdMn2eu+5HHgNmAV8RfPU5t0WQN76xhv0pCfHUWy1oWlQYq0hOyOwCfemqHVIaK/XnMtm4/SHa3znPa++Kmy3FzAkphI3+TY2//dLjtgENSXx1H+ZgZy+6Ywa15tzR/ZEp3fz0lPrfSJUS7S64CoUXYi2tPVNcQA4RwiRAVQAk4G/tFRxIEK0ArhECLEJ0AG3CiHuxTMGuEoI8Rwe8dADv/F6wj0GvO71srAAN7dQ/t+9eTcA9rPkbURmajzF3gWTlrLqdgkReMSoPe7aAIVfrsNptQIQl2XGfPGEdpXXUeQeL2Hr+m85sLsMt7v+3JDRqGfY6F6MndiP7r3qetV2u73ZXU6j0QVXoehCtLqtb6oQKaXDe9/H3ryvSilPtVTxWYVISukG7miQfNDv+ivAKw3uyQcua6HMfn7HVcCcs9nRHOa0eA6f9PT6isIg+KnmdpO38n3feY/ZV6IzhM+6GU3TOHbIwobPjnD8aOMtIxL0FQyOO8jIkRn0nXN5k2WkpJmaFKNodcFVKLoCbWnr/a691uD8feD9pvI2RcTPLodbdIWyvfuoPuXZCdGQkEC3S77X6TY0FXp+2KheHNx7hg2fHeF0buP5n74DMxgxwEXyN6+h12m4DkHlwbEkDhnfKO+UmYNYs+KAWjGuUCiCQsQLkbneTq2hF6KCTz/3HWdNnYwxoX1Dha2lqVheq97axdrVBxrFfNPrdQwf3YtxUwbQvadnoW2B/X9U7PkSgMKP/oGp91AMifUdXs4b1ROj0ahWjCsUiqAQ8UKU6RfmJ9RDc86KSoo2b/Gdh6I31FQsL5dLqydCRqOeUeP6MGHqANIazKllzvwx1d/uxlVRjLvKimXNK2T/4L5GzhaRtGJcoVCENxEvROYwGporXP8Vbu+ivMQB/TtkB9az0ZLDQJzJyIUX92PcpP4kJsc1mcdgSiTryp9xZvljAFQe3Ezl/o0kndf5EcMVCkXXIOJjr5jr9YhCK0T5n37mO+42o3N7Q5pbY++OU81uGBcXb+Seh7/H9MuHNCtCtSQMHEXy+TN855aPX8FZ0eIyAIVCoWgzES9EGSl1PaLichsuV+DbTQeTimPHfFG2dTExZE2Z1Cn1aprGof35LHpmPf9duhO3u3EUBKNRz+XfH4YpPvBV7JkzbsGYYgbAXV1B0dpXg2azQqFQ+BPxQ3OxMQZSk2Ipq7DjdmuUVtjqReXuLPLX1vWGMieMx5jUtk36mvJ4a24u5rsjFj7/UJJ7vKReekyMAb1Rh63a2WZHAn1cAubL76gbotu/kcrzJhHTb2SbnkuhUCiaI+KFCDyLWssqPHMzltLqThcil81G4bq6gBBtdVIIdPdKS0EFn6zaz5EDBfXuj4k1MG5yfy6aOrBVvZ/mSBg4iqThU6jYs85T75pFdLv1zwFvF6FQKBSBEBVCZE6N59gpz9oYS1kNnb2apXjLNlyVnsjUpu7dSB3Wtk3mzrZ7ZU21g/VrD7Ptq2/rDcEZDHrGXNSHid87h6SzzP+0lswZt1J1dKdne/HyYsrWLSN95m1BrUOhUHRtokKIMv32JQqFw4K/k0L296aj07dt6q2l3St3bDnB5x8dpKrCL1S+DkZekMOUmYMbuWEHC0NCMuaZt1Hw3rMAVO76lIShFxE7UA3RKRSK4BDxzgrQcFFr564lqjlzhrLd3tiAej3Z35vW5rKaC5GjN+hY/fbueiLUu38GP/n5JK6+8fwOE6FaEs+9mIRBY3znJR+/gtsZOXvHKBSK8CY6hCiEPaL8z77wHaePOp+4NuxdVEtTu1cC9fYDSkkzce0PRzP/zgn0yGlxi4+godPpMM9aiC7WI5TOktOUftXsHlcKhULRKqJCiDJDFOZHc7ko+KwupE97IykMH5PDFdc27WZtNOqZPPMc7vzVNM4b1bPTt5UwpmSSOf2HvvPSLSux5X/XqTYoFIroJCrmiPwXtXbm0FzJzv9hLyoGICY1hfQLxpzljpbJz7OybeN31FTX34b43JE9mHHl0A4fgjsbyaNnYt37Ffbcg+B2cerVX4HbiTHFTPq0uSQPa9/2GQqFomsSFUKU6b+otawGt1trNsJAMCnwc1LImjYVfRu3PXY4XKxfe5jNXxyt5w2X1T2Zy38wjL4D2z7cF0x0Oj3ply4g/9/3g9sFbs/mjU6rBcsHLwMoMVIoFK0mKoTIFGckMT6GymoHTpcba6WdtCC7MTfEXlpG8ba67d/9h+VatSj1aBGr/7ObYkulL81g0DPpknO4eNpADO3Y+rwjiMnogS4mDs1WVS9dc9oo+WKpEiKFQtFqokKIwBP8tNI7pGUpq+5wISr8ch2ay7PmJ3mIICHHIzSBLkqtrrLz6eoD7Nx6sl65fQZkcOV1IzB3a1tkhs6goQjV4rQ23mhPoVAozkbUCFFmWjzHz5QDHs+5QTlpHVaXpmn1Qvr494bOtigV4MiBQtas2E9Fed3WDHEmIzOuHMrocX3QdcKwYnswpJhxWS2N0o0p4TGEqFAoIouoEaLOXEtULg9RnZsLgN5kwnzxRb5rLS1KddhdfPzefnZuza13TQzrxqwfDCMlBDHy2kLKpBso/fifaM76G+2ljL86RBYpFIpIJoqEyG8tUQe7cOev/bSu3okXY4ivE5DU9PgmxSgpOY5Fz6ynqLCyXtqsHwxj6IgeHWpvsEk8dyJGo5Hiz5fgKq8bjqv5djfaBbM63bVcoVBENuE1E94O/HdqtXTgolZnVTWWDZt85w3XDjW1KFWv11FZYasnQkOGd+eO+6dEnAjVkjxsMn3vXkSPeb/3pVUd/prK/RtDaJVCoYhEokaI/IfmijpwaK5o02bcNZ7y43vnkCwG17s+fEwOV84Z7gvXYzDqcbs1NK9XdmysgcuvPY85t4whITG2w+zsLOL7nEfK6Et955ZP/oWrsiyEFikUikgjaoTIP/BpR/aILF9t8B13+970Joehho3uxbRZgjiTEZezbqO+nL7p3Hr3BEZc0Cuqhq8ypv8QQ+0melVWLJ/8K8QWKRSKSCKK5ojqOytomhb0xt5htVJaG+AUz/xQQ+w2Jx+8s4c9O0750nR6HZMvOYdJ3xuE0+UMqk3hgD4ugawmNtFLHHxhiC1TKELLV8e3sWz3SoqqislMyOCmEVczqe/YUJsVdkSNECWYjMTHGai2ubA7XFRUO0hOCO7QV9HmLeD29HCShwjissz1rlsKKnj7tW8ozK/wpaVnJvD9uaPI6ZvuSajv2R01eDbRm0rFni8BsHz0D0x9zsVgSgytYQpFiPjq+Db+8fVS7C7vpp1Vxfzj66UASowaEDVDczqdrn7w0w4YnvN3UjBPvKjetQO7T/PPv26oJ0Lnj+3Nwvsm14lQlJN5yXwMiZ71W66KEoo/fT3EFikUoWPZ7pU+EarF7rKzbPfKEFkUvkSNEEHHOizYS0sp27vPc6LTkXmRR4jcLjdr39/P269vx27zDLsZjXquvnEkV90wkti4qOl0nhVDfDKZl93uOy/f9RlV3+4KoUUKRegoqipuVXpXJqqEyN9hoTDIPaKiTXXDcinnDiUuM4MKaw2L/7GFzV8e8+VLz0zgx3dfzMgLewe1/kghacgEEoeM951bPngZt73zd81VKEJNZkJGq9K7MlElRNnpddsknPELIhoMLBvq1seYL76IE98Ws+jZrzh+tO7t5pxzs7n95xPp3qtzNqwLVzIvvR19vCdWnrOsgOIv3wyxRQpF53PTiKuJNdSfp441xHLTCBWBpCFRNW6Uk10XKDS3oKKFnK3DVlSMdf8BADS9nm9j+vL5S5vrtmzQwdRLBZO+Nyjs48R1BsakdDIvuZXCVc8DYP36I5KGXoyp95AQW6ZQdB61Dgkd4TWnaRpVjmrKbOVYa8q9Pysot1dgtVVQ7vtUYrV7ju/pN6/d9XYUUSZEyb7j3ILyoJVbtGkzaBpudBw753KOf1w3FBefEMMPfjiagSIraPVFA0nDplCxbwPVR3cCGoWrX6TX7X9BH9OxUdEVinBiUt+xAQuPpmlU2qsorbFSWlPm/en9VFvriU6ZrRyXO3pccKNKiHpm1bkKnymuwuF0ExOE/XwsGzbi0Meyp/s0Slx1gtOzdyrX/WhMyHdODUd0Oh1Zl9/ByX/8HM1ejaM4j+LP3sB82U9CbZpC0alomka1s4bi6lJKqssoriqtO64upaS6lGKv8ESTuLSGqBIiU6yRrPR4Ckuqcbs1zhRV0rtb8tlvbAFboYWCI6fY1etyKuPq3LBHjOnFldePwGg0tHB318aYYsY888cUrn4RAOv2NSQMvpCEAeeH2DKFInjYnXYs1SUUVXk+lqpiLH7HRVUl1DSIVB8M4o0mUkzJpMYlk2JKJiU2kRRTMsmxSSTHJZISl0Ry7Sc2kYN7DgTdhmARVUIEkJOVRKE3+nVuQXm7hWj/mk18nXMlDmOda/jUywYzacY5URWmp6NIGjGNykPbqDr0NQCF779IzoJnMMS373tRKDoLu8uBpbKIgspiCiuLKKwq4oy1gMKqYizVJZTVWINWV3yMiTRTCmmmVO9PzyfV99MrPHFJxBojP1ZlLdEnRN2S2XmoEGjssFCwbj0nFi/FZikizpxJn3lzyZ7S/NbW+3fl8cH/NNxeEdLr4OqbRzF8dK+Oe4AowzNE91NyTx3CVVmGq6IYy5pX6Pb9e0NtmkLho9JeRX5FIWcqLJypKCC/wsKZikLyKwopri5td/kxhhgy4tPIiE8lPT7Nd5zhPU6PTyXNlEpcFIlLazirEAkh9MBLwEjABtwupTzid/0nwELACTwmpVwthDADbwLxQB5wq5Syqpm8GcAhYK+3yBVSyr+19YGa85wrWLeeoy++jNvm6SLbCi0cffFlgHpitGd7Lp99eBBrae2CWM/QW4yrhhtuH8+AYUqEWoshMRXzrDvIf+cpwBOLrmLwWJLOmxhiyxRdCbvLwZnyAvLK832f09Z8zlQUUm5v+3IPvU5PZnwamQnp3k8G5oR0zH7HSbGJYT+CEuS2/jngYqDWa+xqKWWzYfkD6RFdA5iklBOEEOOBp4GrvYZ1B+4GLgBMwAYhxFrgUeBNKeVrQogHgYVCiGXN5B0NLJNS3hWALWelvhDVec6dWLzUJ0K1uG02Tixe6hOiPdtzef8/u3H6RcwGSLCXMTkjjwHDrguGiV2SRDGWpBHTqdj9OQCWNYsw9R6qthdXBJ1KexUny06Ta80j13qGPOsZ8srzKawsRkNrdXk6nQ5zfDrZSWayEjLJSswgPS6NrMQMeqV2Jz0+FYM+KuaKg9LWA8/iadcvlVJaAqk4ECGaCKwBkFJuEUJc4HdtLLBRSmkDbEKII8AI7z1PePN85D0+2kzeMcBoIcQ6oAC4W0p5+mxG2e32JtO7pdW5B+fmV2Cz2dDpdNgsTf8+bJYiX1mffnCgkQilVZ9mxOnPyZn942brbA0Oh6PdZYSS9tifMnUu1cf34CorxF1TSf77L2C+7sFOf1Psyt9BuBCMZ6h21HCq/Ay51tPkWs/4jktqWr8fVow+huzETLolmslONNMtyUy3RDPdErPITEjH2EBoau2PiYnB5XThio5oxkFp64UQfwPOARYJIboB/5JSvtpSxYEIUQrg/826hBBGKaWziWvlQGqD9KbS/NMPAtullJ8KIeYCzwONuh5CiAXAAoAlS5Y0a2xacpwvCneVzUlphZ305DhiM83YmxCj2EzPG7mtxkl5Wf0eU3frEYYWeCIqpF2otjRoL/q4BDJm/ZTC5X8ANGzf7aZy51qSRs8MtWmKMEbTNIqqSzhRlsfxslOcKDvF8bJTFFQG9LLtQ4cOc0IGPZKz6ZHk/SRn0z0pm3RTCnpdVAWaaYTVajUKIb7xS1okpVzkdx6stj4RTzv+DJ65jS+EEN9IKXc3Z1sgQmQF/F2c9F7DmrqWDJT6pVc3kdYw71agypu2Avg9TeD9hS0C2LFjhxYb2/ykXk52ModPeiYYC0psdMtMpu+P5tabIwLQx8XR90dzcdph+b+21ysjtfoM5xZ8hQ4oSetHQnpas/W1hZbsjwTaan/swJHYx8+mEf3ukQAAG6VJREFUbMsqAMrWLSXpnFHEZnb+3FtX/Q7CiYbP4NbcnCkv4GjxCY6WHOd4aS7fleZSaa9qpoTGGPVGeiV3o1dqD3qn9KBXSnd6Jneje3I2sYaYDrU/nElJSXFKKS9oIUuw2voq4G9SyioAIcTneOad2iVEG4HZwH+844Z7/K5tAx4XQpiAOGAoHqeDjcDlwGvALOCrFvK+DrwL/Af4HlBfEdpATnaST4hyC8oZPsjsmwdq6DUXf/5YXn9pU73tGwBE4RZqB4x6TpvUXpMUfqRPuYmqo//DUXgCzWmncNXz9LzlcXTRMc6uCBBN08ivKORo8QmOlRznaPFxjpWcoNoRWOR8g05Pz+Ru5KT2pHdqD3JSetAntSfdkrKiZc6mswlWWz8YWC6EGI0nnulEPO18swQiRCuAS4QQmwAdcKsQ4l7giJRyldc74itvhb+RUtYIIR4DXvd6WViAm6WUlc3kfRB4VQjxM6ASuL2xCa2jfqifOoHJnjK5nodcSVEVr724iZIi79uWDs6/MIfTu4+QbC/xpBljGHH9pe01SeGH3hhL9tX3cOrVX4HbiS3vMKUb3yV90vWhNk3RgVQ7ajhc9C2Hir7lYOERjhUfp8IRWE8nMSaevmk59EvL8fxM701OSndigtzD6eIEs61fCmwBHMAbUsp9LVWs07TWe5GEmh07dmijR49u9vqm3Xk8+bpnAeVokc3vFkxolKcwv5wlL2+h3OoZqtPrdVxz0/kMG92LE8v/w8llbwGQOWEcQx58IGi21zo8RFKX3p9g2l+66b8Uf+HZsRKdnp4/egxTjmh3uWdDfQcdT21vR1qOcajoGIcsxzhhzSOQ9iY1LpmBGX0ZkNGXAem96ZuWgzkhI6zcnyPhO2jI9u3bq8aMGROWWyZH3YJWaN6Fu5bTuWUsXbSVqkrPPyaDUc+cW8Yw+NxuQP0tHzIvvriDre26pI6/msrD27HlHgTNTf6KZ8i57S8YElTUhUjDrbnJLTvN/sLD7C84zIHCw5TZzh54ODE2gUEZfRmQ3peBGZ5PRnxaWImOouOJSiHqYU5ErwO35tkgr8buxBTredQTx4pZ9q9t2Go8c3AxsQZuvO1C+g8yA1B5/ATVJ3MBjzNDxoVjQvMQXQCd3kD2Nfdw6p/3466pwGW1ULDqObrf8BC6KPdginTcmpsTpXnsLzzkE56zLQrV6XT0Te3F4MwBDEjrw6CMfvRO76lERxGdQhRjNNAtM5HTlko0DU5bKunfM5UT3xaz9JWtOOwen39TfAw3/2QsOX3rgpn694bSLxiDwWRqVL4ieMSkZpM1+//If/uPAFQf3UHZ5pWkXfT9EFumaEhBhYXd+QfYdeYAewvkWT3ZEmMTGJzZn8GZAxDmAQzM6Ed8jOfvqXZoS4mQAqJUiMAzPHfau0trbn4FsU6NN1/Z5hOhxOQ4frhgHN16ptS7r2jTFt+xeeJFnWdwFyZx8IWkjr/K59Jd/OWbmHoPwdR7aIgt69pUOarZV3CIXWf2s/vMAc5UFLaYPzkuiXOzzuHc/9/emUdHcV15+OtVC0JSC6EFBMJgeKzCWGZxbBYn8STYHjvxJN7izNiZTJbJyXIyyYyTeJxxjpPZHM9kJokHJx47xoQQ20y8YGJCcNhihCywDYY8NiMQAgktSEi91jJ/VLXUQt1CILVa3X7fOX266tWrqvu6uupX9eq+e8dPZ07JDCoKyjN+bI5ieMhgIRpL7YEmAI4ebeGNF/YTDlndcWPyvPzVF6+luDSvzzqB02cINNjdcl4vvurEDhGK4aVoxacINkhCDdJ6X7T+MSo++yiuMe/vtOsjiWma1J9roK5xH2+dOcDh1vcwTCNh/YKsscwqmc6c8TOYXTKdivxy9YSjuCwyWIgskckGTtSeAjt0T06uh3u/sKSfCAG0v9k76Lhg/jxcWSqb6EjhcLkp/fjXafj5NzAC59G72mh+6UeU3fWgel+URMJ6hHebJXWn9lF3eh+t/vaEdb0uD7PHT6eqbDZVpTOZVKDe7yiGh4wWoixA4OwRoaxsN/d+fgml5flx12nb3StERdcMNABZkQzc+cWU3PoVzqz7PgCBY29zbud6fNerYLPDyblAB3tO7+fNxn3sO3OQkJ44huIVhZOoKpvF/LJZzCieNiyRCbbX7+aXb/+G1kA7xblF3F1126DTaSsyk4wVonyvG4ETrx0fweN1cc/fLKa8In5Xj+b30/nugZ55n/KWSwm5V15N4Qdu59wf1wPQvnUt7VvX4s4vxnfDpxg7N3H+KEViWv3t1DTsZdfJPciWYwmjUOd6criqfA7V5fOoKptJQXb8m7bLZXv9blbVriFsi1+Lv41VtdZYMiVG718yUojOdwT5v1/UkWWLkIHJR++oYtIUH/vqGtiyUdLRHqDAl8MHVwrmVVdwbu9bmLrtyDD1CrLGqfQEqcK3/C66ZQ2R1lM9ZVpnCy0brPxRSowGR4u/jV0n91Jzcg+y9VjCeuVjS6gun0f1xCpE8bR+kaaHk7XvvNgjQlHCepi177yohOh9TMYJUXdXiNWrdvWE7TEwOYyBmethX10Drzy3j0jEEpyO9gCvPGeFU8qq7Q1xV7RQdculEofThRHq7xpsaiHaX1+jhGgAWvxtvHFiD7tO1nG47XjcOg6Hg1nFV1I9oYrqifOYMLZ0xOxr9bddUrni/UFGCVHAH+bZVTW0RAOYOuCoadCJFWHhva3He0QoSiSis+XVgyw5tKenzKeEKOXoXfHTM2udrSNsyejHHwlQc3Iv2+prONB8OG63m9PhZE7JDK6ddDULJ84f9i63wTIut4iWOKIzLrcoBdYoRgsZI0SapvPrp9+kqbETAIcDKhZMoHbPScAKftrRHoi/8ukTaJ3Weh5fIXnTpo6IzYrEuPPHoXX2zzfjzB6VobJGHN3Q2dO4n231NdSeepuI3j/RnMvhZG7pTJZULGBhxVXkZ/X3FB1p7q66rc87IgCvy8vdVbel0CpFqskIITINk5d+9Tb1R3vvtG69cz6hMR7W20J0qrmLEl9OXDGaqPcmhPVVV+NwKnfhVOO74VO0bPgfTO3C9O5+AvXvklM5J0WWpZb32k+y5ehO3mioozPU1W+5w+FgXslMrpt8DQsnzicva3QJd/Q9kPKaU8SSEUK05beS/Xsbe+Y/dPNM5i+cRGNL74na0NzFXbfM6fOOCMDjcTHROEN02F7RItUtNxqIvgdqf32N9WTkdIGhW4NdX/g3Jt73z3iKJqTYypEhEAmy80Qtm4/u4Fj7ibh1KgsmsnTKYq6vXEhRzvAmcRxullYuYnH5VUB6Ra9WJI+0F6K6N+rZ+fsjPfPV11bygRumAVDqy8XtcqDpJm2dQabNKeUW6OM1t+LaYrr+y3oicng8FM6vSkUzFHEYO3dZjyBpHWc59dQD6N3nMAJdnFn3z0y47we4cjI3Uvextnp+d3QHO07UErrgyRDAl13A9ZULWTZlMZWFFSmwUKEYHtJaiA4fbObV9ft75qfPKmHlx+f0jPZ2uZyUF+dxsskKR9/Q3MW86grmVfeetI2vvEr0uamwaq4KcjpKcReMp/STD3D62YcwtTCRtkaaXniU8rsfxJFBydECkSA76mvZfGw777Wf7Lfc43SzcOJ8lk5exIIJc3GqbmRFBpC2QnS6oYPnn6nDNCwPofKKApZdEWDP57/YJxV4RUlfIZox2ddnO+21vdEUfCqawqgme+J0xt/6ZZrX/xCAYP1+Wjb+jOKbv5j2oWZOn2/m1UNb2Hp8F8E4Tz8T88v48NTrWT5lCV6HJbxKhBSZQloKkWHA2id7I2kX+HL48GydmtWvcST/BoJTx5CtdXPlUxuZuaCKN8gF+ifJ0/wBOvb3ZrBVuYdGP3mzPkBkeSPtW9cCcP7t3+MpnkjhkvTzujJNk3ebJRsObWFP4/5+btcep5slk67mxmlLEcXTesQ2mkJBocgU0lKIujrDdNkpvrOy3dz92UW88fBjHPAtxnBaTQp68jjgW8zM/Xsh34qi3dDc18uo4+23MTUrInfulEqyxo8fwVYoLpfC6/6CSNtpuvb9AYC236/G4ytjjFicWsMGSUSPsPPEm2w4tIX6cw39lk/ML+PGaUtZVrl41Hm9KRTJIC2FSNesO0eny8Ed919DSdlYpEf0iFAUw+nmWE5vTptTZ/sKUZ8gp2oQa9rgcDgYf9MX0M41ETx5EDBpfvFHTPj0I2SVj94xYJ3B82w6uo3XjmyjI9jZb/mC8rncPOODzCudmfZdjQrFpZCWQhTl1jvn96T4Drrj3zla5ZZzduPZbnTdwOVyYhoG7XUqrE+64nB7KP3E33PqqQfQzjVhRkKcWfd9yu99GG/x6PIga/Of46U/bWLzsR2ELxh4muXysnzKEm6acQMT8stSZKFCkVrSVohWfFRQFeP9lpfrpCvQP7RJXq6TIo+Hts4Qmm7Q1O5nQnEeXYePEOmwoykU5JM3/coRs10xPLhy8ym789s0Pv0tjJAfvfscp5/97qgRo6aus7x4cBN/OL4LzdD6LCvKKeSj01fw4anXD7r7bXv9bta+8yIt/jbG5fi4Z/7H1EBQRUaQlm43WZqfnHWP0rx1W0/ZjR+fj/uCoMFul1VeUdI71iT6nqgt1ltORVNIW7zFFZTd+R0cHsvtPipG4Zb+715GiobO0/x419N89dV/YvOxHX1E6IrCSXz12s/w41se4WOzPnJJIrSqdk1PnLbWQDuratewvX53UtqgUIwkaflElKUHCJ9t4ehPrLQAJcuX9YwNipfiYfvxVt45YsUta2jqYtHsvkKkoimkN9mTZlJ+94OcXvsIZiTY+2T06e/hHTdxxOw43t7A+oMbqTm5t58H3IxxU7l99koWlM+5rPc/Kn2CIpNJSyGKYoRCnFi9hpLl1uj7CwerRommDQfLhTt09iz+4/UAONxuCubPHxmDFUkje9Isyu76Dmd+9f2+YnTvw0kXoxPnTvHr/a+w+9Rb/ZbNLRHcPnslc0pmDMkBQaVPUGQyaS1EAKGWi6cFiO2aO3W2i7aY3EMFc+fgzs1Jim2KkSVn8mzK7vq2LUYh9K52W4y+h3fc8MelO9V5huf2v8IbJ/f0ewK6unwut89eyYzi4fHiU+kTFJlMWgpRl9fH9mn3ka11IyLyovUrxsc+EXXR3hbjLae65TKKnMlzKLvzO5xZ11eMJnz64WELknrmfDPPvbuBHSdqMc2+ArSo4ir+YvZNXOGbNCz7iqLSJygymbQUIsPhBIeDoCePd7OrmVrXELdLLkpxYQ5ej4twRCdwvptzJ3vj06mwPplHTuUcyu78NmfW/cAWozYan/0uE+79Hp6i8svebnNXC88feJVtx2swTKPPsmsmVPHJubcMuwBFib4HUl5zikwkLYUoFk23HBQGEiKn00HF+DyONXYwxX8aM2KN5citnEx2aclImaoYQXIq51J2x7csMdLC6OfbaFz9EKWf+CbO8VMGtY3t9bt78uZ4XV40I4JxwRPQgvI5fHLOLVw5bnDbHApLKxextHJRT4gflUJBkSmkvRABiTOvxlBRYgnRlf5et17fNSq23GgmOm6m1d/GuMtIoJYzZV7vk5EWtp6MnnmQguX3kFe9csB1Xzu8laf3PoduWvEML/RYm1c6kzvm3oIonnbpDVMoFH3ICCEq8F3c2aCiJA+nqTOt+1RPmYqmMHqJjpuJCkCLv41VtWsALkuMmp7/N4yQHwydjtdXEzp5gNJbv9wnn5FhGLzTdJDNx3awu6G/BxyA2+nmweVfZnbJjCG0TqFQxJL2QuTxuPjgSnHRehV5Tu46tZk83Xp6cufnM3bG9GSbp7hMhnPcTM6UeUz87KM0r3+M0GkriWLwSB0NP/8Gjpv+mgaPkyOtx6lp2MvZi7hD64amREihGGbSUoicLms8Ruyg1YEInGrE8/R/MjnY1FM24dZbcLhciVdSDMhQu80uxnCPm/EUlpJzxzeRW5/h0Im3aMhy05AN/jefuqTtKHdphWL4SUshKvB5eeiHtwyqbse+/fzpX/4dvas38vbW4qv5xseU2+vlsvPkmzy5d92Qu80G4lLHzZimSVe4mxZ/Oy3+Ns52t/ZMt3Rb8x0hOx9VUeKwOmO9Y1g2ZQmF2fk89+4G5S6tUIwAaSlEg6Vp8xaOPr6qJ+eQ5nTzcsl1yLxKzrT6mVyWn2ILk0Oyn1Z+/e4rSQ03Y5omn5hzE0/WrSNi9EardjlczCieyuq3XqAz1GV/znM+1EVH8Dwh/dITxmXrBpNCGhODEaY4c1hYdT2+qdfhLhhPUW5hUn9HhUJhkZZCFNLD/O3L30l4YTANg/rVazi1/jc9ZR5fIb+btgJ5zhKfh159gr/+0LJ+6w/HRXwo2xjq/ofrJf9AtAba45a3+Ns403WWsBYmpIcJaSFCeoSQFiIQCRLQQgS1IMHY+UiQgBbEHwngjwTxh/10RwL9xukA6KbOH0+8GWfPgyPLncVU3ySmFFRwReFkRHElnprf0lm7wa7Rzfktqzm/ZTVZ5Vcyw1fKP5xoQu86hzvfwHdF8LL3PVyc37+NttfXoHe24M4vxnfDpxg7d1mqzVIohsRFhUgI4QR+CswHQsBnpZRHYpb/DfB5QAMekVK+IoQoBn4J5ACNwP1SSv+l1L2YXfEusE+8vpGt9a9xY20z0xtCPXVzp1TSdd9NvLP5IGAJUUeHyeM1q/usv71+N4/XrEYztZ59XFjnYjzx+kY2n9kATj2hnYnYdryGx2tWo9O77k93PUMgEmRxxVXopoFhGhiG9a2bBrqhY5gGmv39zN7n4z6tPLH7V+iGjmZoaIZuf6LTESK6RsTQ0HSrLGJEeubDeoSIHiGsRwhpYZwOZ1yhAPjKhocG9Tslg2x3FuNziygeU8S43CJrOreI8WOs76KcQpxOZ99xOH/2GXIq53L2lZ9gBHu7b0Onj/Q4NgBonS00vfQTTMMgv2rFSDcNsESo6eWf4rSfErXOFppe/ilA2ojRH+pOUrfxJZYbNfic3eg5Pso/8pdpY/+lkG5tTfa1fqB9Oy4MURLHuNuBW6WU9wkhlgDfklLeZi8rA34HXANkAzvs6X8H9kgpnxZCPGA3au1g60op/2Mgm3Zv2mQeXm+9ZHbiYKx3DMFwmIgRIr/boLBL76nbUDGW+o8vYG/bESJBN2bADvfj1HG4NBwOB3neXAC6wt12zDATYuJTOhyQ68m2o4nZy6FPXRMTXTfQzAg4or9pdDvWt9vp6l0Hs2c63j7f9xguTM2NqbtBd2NqHtC85GflkefNw2Vm4zazcJvZ9nQ2Tjw4BvEjGoYlos6Y1B9eI8jk4CEmBQ9RHq7HSXyhBdBxojs8aA53z7fm8KA73Bg4iT2Q5oV/JGDgMw4S/RH07o4EaztwjSm46FZTTTiiEwmFyHKEL2ihA2dWDg53+gzQjV43EwWyHY1tHXPDnf7q6uqEL0iTfa2XUoYu3GeUwXTNXQ/8FkBKuUsIETv4ZhGw095BSAhxBKiy1/mBXWejPX30EuoOKESOUJgpp2Pv+uO3b6/IYfuCbMzWPwHgzNIgq3/3SjdWgjyHO7EWBOgeyCQL58BaoqMNsDR9MU3AcFniYbhAj5k2XJj2PLrbno4KjAvTsL91jy04btA9YMbPD9Vif2L2DgTsz9DYSSlQSo5jEXM9Dcz31jPT04jH0VeUXBi4zBBeM+F5lRwGOltDLQMsHEV4EpTr9ieTGGVtjd+h3odkX+trE+14MEKUD3TEzOtCCLeUUouz7DxQcEF5vLKL1e2HEOJzwOcA1qxZg/PBbwLgwGF6nO5A2IjkXrhOtf3xOj3+sKHlJrqb9DrdfoDB1BmIsGbk4jSIt42onQOuP8T9Axim4dZMvf/tlmHdqXvdzkFtZyDCmtHvt44yHNsfqX0MhKlNznVwLV0XluPA4fYkff/xbYrkOuL8P1Jp06UQ1oxc70A3Y27vqG/DYBmNbQ2HwzlCiNiXrE9IKZ+ImU/2tT4hgxGiTmBszLzTNizesrHAuZjyQJyywdTth/2DPRFvmUKhUCgujpQDZitI9rU+IYPJj70TuAnA7jfcF7NsN7BUCJEthCgAZgH7Y9cBVgLbL7GuQqFQKEaWZF/rEzIYZ4WoJ0UV1iuQ++0dH5FSvmR7R3wOS9R+IKV8QQhRCvwCSwlbgHuklN2XUndAoxQKhUIxrCT7Wj/Qvi8qRAqFQqFQJJPBdM0pFAqFQpE0lBApFAqFIqWkTYifi436TYE9e+l1UXwPWAX8CGsk8SYp5cOJbLZfBF523SHavRj4VynlCiHElcDTWD7j+4EvSSkNIcR3gZvtfX5NSrk7WXWHoQ1XAy8Dh+3Fj0sp143GNgghPMD/AlOALOAR4EAybErWMUjQhgbS5xi4gJ8BAms0z/1Y70OG3aZknweZRDo9EX0MyJZSXgs8APwwVYYIIbIBpJQr7M/9wP8A92AN8FpsXyAT2TzUupdr998DP8ca7QzwGPCglHIp1sl4m7395cBi4C7gJ0muO9Q2XA08FnMs1o3iNtwLtNrbXAn8OIk2JesYxGtDOh2DPweQUl4HPGTvI92OQcaRTkLUZ9QvVviIVDEfyBVCbBJCbBFCLAOypJRHpZQm8BrwIeLYLITIH4a6l8tR4PaY+Wpgqz29EfiwbccmKaUppTwBuIUQ45NYdzjacLMQYpsQ4kkhxNhR3IbngH+MmdeSaFOyjkGiNqTFMZBS/gZ7YDxQCTQl0aZkngcZRToJUdxRvymyxQ88CnwE+ALwlF0WJdEIY90u6xxi3cvCdqGMxBQ5bIEbyI5oebLqDrUNu4FvSimXAceA747WNkgpu6SU5+0L9fPAg0m0KSnHIEEb0uYY2G3QhBC/AP7bbkNaHYNMJJ2EaKBRvyPNIeBZ+67mENYfKzZjW6IRxs44ZZdTd7iIDaJ2sVHRyao7VP5PSlkXnQYWDJNdSWmDEGIS8DqwWkr5yyTalLRjEKcNaXUMAKSUfwXMwHpflJMkm0byPEhr0kmIBhr1O9J8BvsdjhBiApALdAshpgkhHFhPStuJY7OUshMID7HucLFXCLHCno6Oit4JfEQI4RRCTMYS/JYk1h0qrwkhovk1PgTUjdY22IP/NgH/IKX8X7s4rY5Bgjak0zH4tBDiW/asH0ss3kynY5CJpI3XHNad1o1CiD/SO+o3VTwJPC2E2IHlEfMZrD/0GsCF1QdcI4SoJb7NXxhK3WFsx98BPxNCeIGDwPNSSl0IsR14A+tG5UtJrjtUvgj8WAgRBs4An5NSdo7SNnwb8AH/KISIvmf5KvBfaXQM4rXh68B/pskxWA88JYTYhhUf+2v2ttP9PEhrVGQFhUKhUKSUdOqaUygUCkUGooRIoVAoFClFCZFCoVAoUooSIoVCoVCkFCVECoVCoUgpSogUCoVCkVKUECkUCoUipSghUigUCkVK+X9GG6DszxKzXgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax1 = plt.subplots()\n",
    "ax2 = ax1.twinx()\n",
    "ax1.plot(calibrator.lab.time, calibrator.lab.NO3.concentration[0], c=sns.color_palette(\"deep\", 10)[0], lw=3, label='NO3')\n",
    "ax1.scatter(NO3_t, NO3, c=sns.color_palette(\"deep\", 10)[0], lw=1)\n",
    "ax1.plot(calibrator.lab.time, calibrator.lab.SO4.concentration[0], c=sns.color_palette(\"deep\", 10)[1], lw=3, label='SO4')\n",
    "ax1.scatter(SO4_t, SO4, c=sns.color_palette(\"deep\", 10)[1], lw=1,)\n",
    "ax1.plot(calibrator.lab.time, calibrator.lab.CH4.concentration[0], c=sns.color_palette(\"deep\", 10)[2], lw=3, label='CH4')\n",
    "ax1.scatter(CH4_t, CH4, c=sns.color_palette(\"deep\", 10)[2], lw=1,)\n",
    "ax2.plot(calibrator.lab.time, calibrator.lab.Fe2.concentration[0], c=sns.color_palette(\"deep\", 10)[3], lw=3, label='Fe2')\n",
    "ax2.scatter(Fe2_t, Fe2, c=sns.color_palette(\"deep\", 10)[3], lw=1,)\n",
    "ax2.plot(calibrator.lab.time, calibrator.lab.CO2.concentration[0], c=sns.color_palette(\"deep\", 10)[4], lw=3, label='CO2')\n",
    "ax2.scatter(CO2_t, CO2, c=sns.color_palette(\"deep\", 10)[4], lw=1,)\n",
    "ax2.grid(False)\n",
    "ax1.grid(lw=0.2)\n",
    "ax1.set_ylim(0, 2.5e-3)\n",
    "ax1.set_xlim(0, 40*24*60*60)\n",
    "ax2.set_ylim(0, 25e-3)\n",
    "ax1.legend(frameon=1, loc=2)\n",
    "ax2.legend(frameon=1, loc=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernel_info": {
   "name": "python3"
  },
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  },
  "nteract": {
   "version": "0.11.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
